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ABSTRACT 


The  FEE,  or  Flat -Echo  Estimator,  represents  an  entirely  new  solu¬ 
tion  to  the  problem  of  spectrum  recovery  from  a  given  finite  set  of  error - 
free  covariance  samples.  In  concept  it  is  based  squarely  on  frequency- 
domain  network-theoretic  ideas  and  its  principal  features  are  the  following. 

1)  Numerical  robustness,  as  defined  in  this  report,  can  be  assigned 
in  advance  by  the  adjustment  of  a  single  parameter. 

2)  It  is  possible  to  "tune"  the  estimator  to  produce  selective  ampli¬ 
fication  of  a  finite  number  of  desired  frequencies  without  impairing  either 
its  interpolatory  character  or  its  numerical  robustness.  Enhanced  reso¬ 
lution  is  to  be  expected. 

3)  For  the  parameter  setting  which  yields  maximum  robustness,  the 
FEE  coincides  with  the  MEE,  the  Maximum -Entropy  Estimator. 

In  addition  to  providing  a  rather  leisurely  account  of  some  of  the 
network  ideas  behind  the  FEE,  this  research  report  also  contains  two 
appendices,  A  and  B.  The  first  supplies  an  explicit  and  compact  formula 
for  the  class  of  all  interpolatory  spectral  densities,  while  the  second  is 
devoted  to  the  design  of  minimum -degree  regular  all-passes  with  prescribed 
phase  at  a  given  finite  set  of  frequencies. 
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I.  INTRODUCTION  AND  PRELIMINARY  RESULTS 

Suppose  that  we  have  a  discrete-time  random  process  where  "t"  can 
traverse  all  positive  and  negative  integers.  Let  us  also  suppose  that  the 
process  is  zero-mean  and  second-order  stationary.  Then, 

^K*t+k^  =  C(k)  .  Ikl  =  O-^oo  ,  (1) 

is  the  associated  covariance  function.  Clearly,  for  all  integers  k, 

C(-k)  =  E(x^x^_j^)  =  E(x^^j^x^)  =  ^(k)  .  (2) 

As  is  well  known  [l]» 

C(k)  =  ^  /"e‘-’*'®dF(0)  (3) 

-TT 

where  F(0)  is  monotone -nondecreasing,  continuous  from  the  left  and  of 
bounded-variation: 


VarFs/  dF(0)  =  2ttC(0)  <  <«  . 
-TT 


Further,  the  derivative 


K(0)  = 


dF(0) 

d0 


>  0 


(4) 


(5) 


exists  for  almost  all  0,  is  L^  over  -tt<  9  <  TT  and  defines  the  spectral  density 
of  the  process.  Evidently,  if  F(0)  is  absolutely  continuous, 

C(k)  =  ^  /’^e’-’*"®K{9)d0  .  (6) 

-TT 


Let  us  assume  first  that  F(0)  is  absolutely  continuous  and  that  K(0) 
satisfies  the  Paley-Wiener  criterion. 


TT 

/  fnK(9)d9>-<»  .  (7) 

-TT 


a  =  complex  conjugate  of  a. 


I-l 


Then,  there  exists  [  1  ]  a  imique  function 


B(z)  =  Y,  *  ^0  ^  °  ' 

r=0 

which  is  analytic  and  devoid  of  zeros  in  )z|  <  1  and  is  such  that^^^ 

r=0 

and<2) 

K(0)  =  |B(ej®)l^  ,  a.e. 

Moreover,  x  can  be  realized  as  the  output  of  a  causal  digital  filte 

^  (4) 

square -summable  impulse  response 

h  =  b  ,  r  =  0->  oo 
r  r 

driven  by  an  appropriate  white  noise  source  sequence  u^: 

00 

S  Vt-r’  (-«<t<*)  . 

r=0 


E(u  )  =  0,  E(u  u  )  =  6  ,  (-<»<r,m<oo) 

'  r  ’  '  r  m  rm  '  '  ’ 

00 

^^^t^t+k^  ■  "  S  ^r+k^r' 

r=0 


7^ 


B(z)  is  of  class  H^. 


(3) 

(4)  00 


id  id 

B(e'*  )  =  limit  B(re''  )  exists  for  almost  all  6. 
r-*  1  -0 


S  l^rl  =  S  l^rl  ' 

r=0  r=0 


(8) 

(9) 

(10) 
with  a 

(11) 

(12) 

(13) 

(14) 
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We  can  now  pose  and  solve  several  problems  in  spectral  estimation. 

Problem  1.  It  is  known  that  F(0)  is  absolutely  continuous  and  that  the 
spectral  density  K(9)  =  dF(0)/d6  satisfies  the  Paley -Wiener  condition  (7), 
Under  the  assumption  that  the  n+1  covariance  samples  C(0),  C(l),  . . . ,  C(n) 
are  known  absolutely  and  without  error,  determine,  if  it  exists,  a  compatible 
K(0)  which  maximizes  h^  in  the  input-output  model  (12). 

Solution.  Let  the  numbers  C(k),  k=n+l-»co,  constitute  any  admissible 
completion  of  the  given  data  C(0),  C(l),  . . . ,  C(n)  and  let  K(0)  be  the  corres¬ 
ponding  spectral  density. 

Let 

TT 

dk  =  ^/  e-j^®inK(0)d0  ,  (15) 

-TT 


k  =  0  ^  00,  let  =  1  and  let  the  coefficients  a^^,  k  =  1  oo,  be  determined  from 
the  equation 


exp(  S  E  • 

r=l  k=0 


(16) 


Then  [l], 


hj.  =  arexp(do/2)  ,  r  =  0-»oo  . 


(17) 


Theorem  1.  To  construct  an  admissible  completion  C(k),  k  =  n+l  -» oo, 
which  maximizes 

is  equivalent  to  choosing  K(9)  to  maximize 

do  =  ^  /  fnK(0)d0  .  (19) 

-TT 


But  [9  ]  this  leads  immediately  to  the  MEE,  the  unique  maximum -entropy 
estimator,  Eqs.  (73)  and  (74). 

Comment  1,  What  does  it  mean  to  maximize  h  and  is  it  also  possible 

o  ^ 

to  select  spectral  density  estimators  which  maximize  other  meaningful 
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functionals  defined  on  the  impulse  response  sequence  {h^}?  These  questions 
deserve  to  be  studied  and  one  possible  approach  is  to  combine  (17)  with  the 
following  explicit  formulas  for  the  and  dj^'s  [2]; 

d2  ^3 

d^,  ^  d^+  d^d^+ (20) 


and  in  general, 


a  =  — r  det 
r  ri 


‘^l 

-1 

0 

...  0 

2^2 

-2 

•  •  •  0 

(r-l)dr.i 

(r-2)dr.2 

(r-3)d^_3 

.  .  .d^,  -(r-1) 

rd 

r 

(r-l)dr.i 

(r-2)dr_2 

...2d2,  d^ 

.  r  >  1. 

(21) 


=  ‘^r+^l‘^r.l  +  *--'  •  (^2) 

Problem  2.  Let  F(6)  and  K(6)  be  subject  to  the  same  constraints  as  in 
problem  I  and  assume  once  again  that  the  covariance  samples  C(0),  C(l),  •••, 
C(n)  are  known  absolutely  and  without  error.  Provide,  if  feasible,  a  geo¬ 
metric  time-domain  interpretation  of  the  process  whereby  every  admissible 
completion  C(k),  k  =  n+l-*oo,  is  generated. 

Solution.  As  is  well  known  [  2  ],  under  the  given  hypotheses, 


C(0) 

C(l) 

.  .  .  C(k) 

C(-l) 

C(0) 

.  .  .  C(k-l) 

C(-k) 

C(-k+l) 

.  .  .  C(0) 

for  any  possible  admissible  extension  C(k),  k  =  n+l  where 

det  A  =  I A I  =  determinant  of  square  matrix  A. 

The  Paley-Wiener  condition  (7)  precludes  the  possibility  0. 
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C(-k)  =  C(k)  .  (24) 

Thus,  at  the  first  step,  C(n+1)  must  be  chosen  from  the  set  of  values  §  which 
satisfy  the  inequality 


C(0) 

C(l) 

.  .  .  C(n) 

C(-1) 

C(0) 

.  .  .  C(n-l) 

C(n) 

•  •  •  « 

C(-n) 

•  •  •  •  • 

C(-n+l) 

.  .  .  C(0) 

•  •  • 

C(l) 

f 

C(-n) 

.  .  .  C(-l) 

C(0) 

>  0  . 


(25) 


However,  it  is  not  difficult  to  show  [  2  ]  that  this  set  is  the  interior  of  a 


circle  ©  with  radius 
n 


n 


n  A 


n-l 


(26) 


and  center 


>  =ilil 

=  n  A 


n-1 


h-1 


C(1) 

C(2) 

C(n) 

0 

C(0) 

C(l) 

.  .  .  C(n-l) 

C(n) 

(27) 

C(-n+l) 

C(-n+2) 

.  .  .  C(0) 

C(l) 

In  fact,  if  Aj^r  denotes  the  minor  of  Aj^  that  is  obtained  by  striking  out 
its  first  column  and  last  row,  it  follows  from  a  fundamental  determinantal 

— Tt)  - 

result'  '  applied  to  that 


(28) 


But  as  is  clear  from  (25)  and  (26), 


(29) 


If  D  is  any  determinant,  and  S  is  the  second  minor  obtained  from  it 
by  striking  out  its  ith  and  kth  rows  and  its  jth  and  f  th  columns,  and  if  we 
denote  by  A^j  the  cofactor  of  the  element  which  stands  in  the  ith  row  and  jth 
column  of  D,"^  then  [3  ], 


A..  A., 

ij  t£ 

A,  .  A,  , 
kj  kf 


=  (.l)i+j+k+i  Ds 


(29a) 


1-5 


1 


and  (Z8)  yields. 


n 


^n-1  “ 


(30) 


since  j  both  positive. 

Once  having  chosen  C{n+1)=  §  consistent  with  (30),  C(n+2)  is  then 

selected  from  the  interior  of  a  circle  of  radius  r  ,  ,  =  A  ,  ,/A  and  center 

n+ 1  n+ 1  n 

given  by  (27)  with  n  replaced  by  n+1,  etcetera.  Evidently,  except  for 
r^  and  5^  which  are  uniquely  determined  by  the  prescribed  data  C(0),  C(l),  . . . , 
C(n),  all  successive  rj^  and  depend  on  the  particular  rule  used  to  pick 
C(k+1)  from  the  interior  of  k  =  n+l  -*ca. 

Theorem  2.  Under  the  conditions  prevailing  in  problem  1,  the  following 
is  true. 

1)  Let  C(k),  k  =  n+l-»oo,  denote  any  admissible  completion  of  the  initial 

covariance  data  C(0),C(1),  •- • ,  C(n)  and  let  K(0)  be  the  corresponding  spectral 

i  8  2 

density  estimator.  Let  K(9)  =  lB(e‘^  )|  where  B(z),  Eq.  (8),  is  the  associate 
Wiener-Hopf  factor.  Then,  the  sequence  of  radii  r^^,  k  =  n-*i»,  generated  by 
this  completion  is  always  monotone-nonincreasing  and 


limit  r.  =  limit  -r - 

k-»«)  k.-»oo  ^k-1 


>  0 


(31) 


In  addition,  for  the  sequence  of  centers,  k=n->oo. 


lim  sup  I  lb^< 
k-*oo 


(32) 


2)  The  maximum -entropy  estimator  is  the  (unique)  spectral 

density  that  is  obtained  by  identifying  each  C(k+1)  in  the  completion  with 

the  center  of  circle  ©,  : 

-  k 

C(k+1)  =  ,  k  =n-»oo  .  (33) 

Under  this  ME  rule, 

limit  ?  ^  =  0 
k-»oe 


(34) 


and  all  circles  ©j^  possess  the  same  maximum  possible  radius  r^: 


r 

n 


(35) 


k  =  n-»oo  , 

Proof  (in  outline).  Given  any  admissible  completion  and  its  corres- 
ponding  spectral  density  estimator  K(6)=  |B(e'^  )|  where  B(z)  is  the  associate 
Wiener -Hopf  factor,  it  is  well  known  [4  ]  that 

^k  2 

limit  - -  =  b  >  0  .  (36) 

k->00  ^k-1  ° 


Furthermore,  from  (28)  with  n  replaced  by  k. 


k+l 


“k+1 

A,. 


k-1 


*k+r 


'S 


=  r, 


(37) 


k-1 


and  it  is  clear  that  the  sequence  r^^,  for  k  =  n-»oo,  is  monotone-nonincreasing. 
Thus, 


=  limit  TT - 

k-»<»  ^k-1 


r 

n 


(38) 


with  equality  iff 

k=n-»oo  .  (39) 

Since  K(9)  is  integrable, 

limit  C(k+1)  =  0  (40) 

k-»<» 


which  implies  that  all  circles  ©j^  contain  the  origin  as  an  interior  point  for 
k  sufficiently  large.  Hence,  for  large  enough  k, 

?kl<>^kl^n 


so  that 


limit  sup 
k->oo 


^kll 


limit 

k-»l30 


-  < 


r  = 


n 


n-1 


(42) 


To  establish  2),  observe  that  (39)  and  (29)  lead  to  the  conclusion  that 
b^  assumes  its  maximum  possible  value 


(b.) 


o'MAX 


(43) 


iff  for  every  k  =  n  -♦<»  , 


§  =  C(k+1)  =  . 


(44) 


the  center  of  circle  ©^.  But  according  to  (11)  and  (18),  b^  =  exp(dQ/2)  so 
that  is  necessarily  paired  with  the  estimator  K(0)  =  Kj^^(6)  which 

maximizes  the  entropy 

do=^/  fnK(e)d0  .  (45) 


Lastly,  since  the  ME  completion  rule  identifies  with  C(k+1)  for 
every  k>n,  it  is  evident  that 

limit  §.  =  limit C(k+1)  =  0  ,  (46) 

k-*oo  k-»oo 

Q.E.D. 

Comment  2.  Theorem  2  suggests  that  the  ME  estimator  Kj^^(6)  achieves 
good  performance  because  its  completion  rule  possesses  a  kind  of  optimal 
robustness  against  loss  of  positive  ^definiteness.  Or,  stated  somewhat 
differently,  since  each  additional  covariance  sample  is  chosen  at  the  center 
of  its  circle  of  admissibility,  it  appears  plausible  to  conjecture,  on  the 
grounds  of  symmetry,  that  the  ME  rule  possesses  a  fair  degree  of  self¬ 
correction  in  the  face  of  (unsystematic)  sample  error.  Admittedly,  these 
considerations  are  heuristic  in  nature  but  they  seem  to  be  reinforced  even 
more  by  the  realization  that  the  initial  data  C(0),  C(l),  . . . ,  C(n)  must  actually 
be  estimated  from  samples  of  the  random  process  x^  itself. 

Problem  3.  Under  the  assvimption  thaf(  C(k)  is  real,  place  the  general 
problem  of  spectral  estimation  from  a  prescribed  set  of  error-free  covariance 
samples  C(0),  C(l),  ...,  C(n),  in  a  network-theoretic  setting  involving  positive- 
real  functions. 
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Solution.  First,  it  is  easily  seen  that  C(k)  real  implies  that  the 


function 


Z(z)  =  C(0)  +  Yt  2C(k)z‘ 


is  positive-real  (pr)  in  |z|<  1,  i.e.,  Z(z)  is  real  for  z  real  and 


ReZ(z)  s  >0  ,  lzl<l  . 


The  proof  is  quite  simple. 

In  fact,  it  follows  readily  from  the  reality  of  C(k)  and  (3)  that^ 

1  r'^  ej®+z  ,  , 

=  ^  •  l^l<  »  • 

-TT  ® 


Hence, 


,  TT  1  -  I  z  1  ^ 

ReZ(z)=-^/  — - -dF(0)>O.  lzl<l.  (50) 

I  ^  I  — 


since  dF{9)>0,  Q.E.D. 

It  can  be  shown  [5  ],  that  for  almost  all  0,  the  "interior”  radial  limit 
Z(e'’®)  s  Umit  Z(rej®)  (51) 


exists  and  that 


r-»l-0 


K(0)  =  ReZ(e''®)  . 


However,  even  more  is  true.  Let  Fg(0)  denote  the  singular  part'^'  of 
F(0)  and  let  the  countable  set  of  discontinuities  of  F(0)  occur  at  0=0^^  with 
jumps  p^>0,  i  =!-♦<».  Then,  (49)  may  be  rewritten  as 

J  0.  TT  . 

zt'l  =  W  ^  ^  .f,  Oi  ^  ^ 

_  e*'  -z  1=1  I  -rr  -z 

-TT  e  -z 


+  Zj,(z)  +  Z^(z) 


.  (54) 


^  ^The  interchange  of  summation  and  integration  is  justified  by  dominated 
Convergence, 

19) 

'  'Fg(0)  is  monotone -nondecreasing  and  its  derivative  equals  zero  almost 
everywhere  [6  ]. 
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V 


Clearly, 


Pi  -  j  Si 

— —  =  limit  (1-ze  )  Z(z)  ,  i=i-><»  , 


(55) 


and  the  pJs  emerge  as  residues  of  Z(z)  at  its  poles  z^=e  ^  (all  simple)  on  |z| 
These  poles  are  revealed  by  the  presence  of  the  Foster  reactance 


"h  II  Pi  ^ 


jBi 

e  +  z 


(56) 


i=l  e  -z 


Of  course, 


ReZp(ej®)  =  0 


(57) 


and  since  F(0)  is  of  bounded-variation, 

00 


Epi<“® 

i=l 


(58) 


so  that  the  infinite  sum  in  (56)  is  well-defined. 

Thus,  given  C(k),  k=0-»n,  the  problem  of  spectral  estimation  is  equiva¬ 
lent  to  that  of  generating  all  pr  functions^ 


Z(z)  =  C(0)  +  2C(l)z  +  . . .  +  2C(n)z’^  +  0(z"‘^^)  . 


(59) 


We  shall  now  present  an  elementary  solution  based  on  the  synthesis  of  ideal 
TEM  transmission-line  cascades  [7,8], 

The  Interpolatory  Cascade.  With  the  given  data  of  real  numbers  C(k), 
k=0'*n,  construct  the  real  symmetric  (n+l)x(n+l)  matrix 


T  = 
n 


C(0) 

C(l) 

.  .  .  C(n) 

C(1) 

C(0) 

.  .  .  C(n-l) 

C(n) 

C(n-l) 

.  .  .  C(0) 

(60) 


TTCT 


Phrased  differently,  find  all  pr  function  Z(z)  such  that  d  Z(0)/dz  , 


k=0'»n,  are  prescribed  real  numbers. 
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An  admissible  completion  C(k),  k  =  n+l-*i»,  exists,  if  and  only  if  T^> 

Two  cases  must  be  distinguished  when  this  criterion  is  satisfied. 

Case  1:  A  =detT  >0.  Then,  rank  T  =n+l  and  all  interpolatory  pr 
n  n  n  — 

solutions  Z(z)  are  constructed  in  the  following  manner. 

Step  1.  Let  Tj^  denote  the  (k+l)x(k+l)  Toeplitz  matrix  obtained  by  sub 
stituting  k  for  n  in  (60),  let 

A^  =  detTj^,  k>0  ,  (61) 

and  let  Aj^/  denote  the  minor  of  Aj^  obtained  by  striking  out  its  first  column 
and  last  row.  Let  N  denote  the  2-port  reactance -cascade  (Fig.  1)  composed 
of  n+1  ideal  TEM  lines  all  of  the  same  1-way  delay  T>0  and  positive  charac¬ 
teristic  impedances, 

Ro  =  C(0)  (62) 


and 


l+(-l) 


k-l  V 


^k"^k-l  • 


*k.l 


k-l 


,  k=l-+n 


'k-l 


(63) 


Then, 


under  the  identification 


(12) 


z  = 


.-2pT 


(64) 


every  admissible  interpolatory  solution  Z(z)  is  obtained  as  the  impedance 
seen  looking  into  the  input  port  of  N  with  its  output  port  closed  on  an  arbitrary 
pr  function  W(z).  (In  Fig.  I,  terminals  cd  are  connected  to  ab.) 

^  ^  ^ ^  1  and  dejiote  the  kxk  i^ntity  and  zero  matrices,  respectively. 

For  any  matrix  A,  A,  a'  and  A*(  =  A^)  denote  its  complex  conjugate,  transpose 
and  adjoint.  If  A(=a‘)  is  kxk  and  hermitean,  A>0j^(>0jj)  means  that  it  is 
positive -definite  (nonnegative -definite).  As  usual  x  =  column-vector. 

(12) 

'  'p=a+juu  is  the  standard  complex-frequency  variable  used  in  network 
theory. 


I-ll 


The  quantities 


R.  “R.  4 

-  _  ^  1  —I  1  \ 

■  R. +R.  , 
k  k-l 


k-1 


^k-1 


k  =  1  -»n 


(65) 


are  the  n  junction  reflection  coefficients.  Their  physical  significance  is  not 
difficult  to  grasp.  Indeed,  let  the  terminated  structure  depicted  in  Fig.  1 
be  excited  by  an  input  current  source  i(t)  =  6(t).  This  impulse  immediately 
launches  a  voltage  impulse  C(0)6(t)  of  strength  C(0)  into  line  0  which 
then  travels  towards  line  1. 


After  T  secs.,  a_  reaches  the  first  junction  separating  line  0  and  line  1 

(13) 

and  part  of  it  is  reflected  and  part  of  it  is  transmitted.'  If  a^^  and  a^^  de¬ 
note  the  strengths  of  these  reflected  and  transmitted  impulses,  respectively, 
conventional  transmission  line  theory  yields. 


(66) 


and 


(67) 


In  the  same  way,  after  T  more  secs.,  a^^  reaches  the  second  junction 
and  is  partly  reflected  and  partly  transmitted  to  generate  a^^.  and  a^^  where 


(68) 


and 


(69) 


In  general  therefore,  upon  reaching  line  k,  a^j^  is  divided  into  an  impulse 
of  strength  reflected  back  into  line  k-l  and  one  of  strength  aj^^  transmitted 
into  line  k.  Clearly, 


^^(k-l)t^ 


^®(k-l)t'^ 


1-Sj^  ,  k=l-»n  . 


(70) 


THT 


It  is  convenient  to  label  an  impulse  by  its  strength. 
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Moreover,  since  interaction  with  the  passive  load  W(z)  commences  only 

after  the  lapse  of  (n+l)T  secs.,  the  input  voltage  response  v(t)  over  the  time 

interval  0<t<  2(n+l)T  is  determined  solely  by  the  {n+l)-line  cascade.  We  can 

2 

now  relate  the  radii  r^^  to  the  squared  amplitudes  . 

For  according  to  Eqs.  (26),  (28)  and  (65), 


1-Sj^  =  l 


/  A,/  - 


'k-l 


'^k-2^k 


‘k-1 


k-l 


k=l-»n 


(71) 


A 


-1 


3  1 


Thus,  since  rQ  =  C(0)  =  a  comparison  of  (71)  with  (70)  immediately  gives 
the  important  physical  connection  formula. 


r 


2 


_^kt 

k  CM 


k  =  0->n 


(72) 


Obviously  then,  since  successive  transmitted  impulses  can  never  in¬ 
crease  in  strength,  the  corresponding  sequence  of  successive  radii  of  the  circles 
of  admissibility  must  be  monotone -nonincreasing. 

It  follows  from  this  interpretation,  that  to  avoid  any  decrease  in  circle 
radius  for  k>n,  it  is  necessary  that  the  last  transmitted  impxilse  in  line  n  s\iffer 
no  reflection  upon  reaching  the  load  W(z).  But  this  is  only  possible  if  W(z)  = 

R^,  the  characteristic  impedance  of  the  last  line  in  the  (n+l)-line  cascade. 

Consequently,  ^^^^,(9)  is  the  real  part  of  the  driving -point  impedance  of 
the  n+l-line  structure  N  closed  on  W(z)=R^.  Explicitly  [9,  10,  Appendix  A], 


where 


*SvIe(®)  = 


1 


P 

n 


(z) 


1 

•JW-l 


C(0) 

C(1) 

.  .  .  C(n) 

C(l) 

C(0) 

.  .  .  C(n-l) 

C(n.l)  C(n-2)  .  .  .  C(l) 


n 


z 


n-1 

z 


1 


(73) 


(74) 
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(14) 

is  a  degree  n  strict -Hurwitz  polynomial. 

Case  2:  A  =detT  =0.  Now  rank  T„<  n+1  and  there  exists  a  li  <  n  such 
n  ti  w 

that  A  =  det  T  >0  but  A^_|_ ^  =  det  =  0.  The  interpolatory  pr  function 

Z(z)  is  iinirino,  is  Foster  and  is  of  degree  |i+l  =  rankT^.  Thus,  according  to 
(53)  and  (54),  Z^(z)  =  Z  (z)  s  0  and 

S  <x 

u+i  je. 

i=l  e  ^-z 

j®i 

The  poles  z^^  =  e  ,  and  the  residues  i=  l-»li  +  l, 

1)  Find  a  real  nontrivial  column-vector  solution  x  =  (xq,x^,  . 
of  the  homogeneous  equation 

T^,+iJ=2  •  (’‘■I 

(There  exists  only  one  such  linearly  independent  x. ) 

2)  The  ti  +  1  zeros  of  the  polynomial 

x(z)  =  Xq+x^z+...+x^^jz‘"+^  (77) 

j®i 

are  all  of  unit  modulus  and  are  precisely  the  desired  poles  z^=  e  ,  i=  1-»U  +  1. 
j^i 

3)  With  the  e  's  known  from  2),  the  quantities  p^,  p^,  •  •  • ,  are 

determined  as  the  unique  solutions  of  the  linear  system  of  the  Vandermonde  type, 

-jke 

C(k)  =  ^  P.e  "  ,  k  =  0-»n  .  (77a) 

w  1 

i=l 

A  much  more  enlightening  physical  interpretation  of  case  2  is  available. 
Namely,  since  0  implies  that  T^>  then 

Aj^>  0  ,  k  =  0-»|J  .  (78) 


(75) 


can  be  found  in  three  steps 


.(15) 


^  0,  Iz]  <  1.  Incidentally,  observe  that  Pn^z)  Is  the  determinant 
of  the  ma'trix  that  results  by  replacing  the  last  row  in  Tn  by  z*^,  zn-1 . .  1. 

This  three-step  rule  goes  back  to  F.  Riesz  [2]  but  its  use  in  the  context 
of  spectral  estimation  is  of  recent  vintage  [l  1  ]. 
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Hence,  from  (71), 


Sj^<  1 


k=  l-*[i  . 


However,  since  s^^j=l;i.e.,  ±1  and  consequently,  in  view 

of  (65),  =  00  or  0.  In  other  words,  the  interpolatory  Z(z)  is  the  input 

impedance  of  a  cascade  of  4+1  TEM  lines  of  characteristic  impedances 
Rq,  R^,  •  •  • ,  R|^  closed  on  either  an  open  or  short  circuit  (but  not  bothi ).  This 
explains  why  the  process  terminates  and  why  Z(z)  is  unique.  Alternatively, 
we  can  say  that  the  choice  of  covariance  sample  C(|a  +  1)  has  been  made  from 
the  boundary  of  the  circle  of  admissibility  and  as  a  result  the  transient 
voltage  impulse  transmitted  into  line  4  +  1  must  be  zero.  This  situation  is 
realized  only  because  the  load  that  terminates  line  4  is  either  an  open  or  a 
short.  (We  conclude  that  the  two  associated  real  sample  values  C(4+l)  are 
the  two  intersections  of  ©^  with  the  real  axis.)^^^^ 

In  the  next  section  we  shall  undertake  a  much  deeper  study  of  the 
properties  of  the  spectral  estimator,  especially  with  regard  to  its  dependence 
on  the  choice  of  pr  load  W(z).  Nevertheless,  for  now  we  can  assert  that 
any  selection  of  covariance  sample  made  from  the  boundary  of  its  circle  of 
admissibility  immediately  forces  Z(z)  to  be  a  unique  Foster  function  of 
finite  degree  and  C(k)  to  be  a  finite  stun  of  pure  tones  (Fg(9)  =  K(0)  =  0): 

,  -jk9. 


'  'Since  the  data  C(k),  k=  0-»4,  are  real,  the  center  of  ©  lies  on  the 
real  axis.  ^ 


n,  "TUNED"  SPECTRAL  ESTIMATORS 

From  this  point  on  it  is  assximed  that  the  given  data  C(k),  k=  0->n,  are 
real  and  error-free,  that  is  positive -definite  and  that  F(0)  is  absolutely- 
continuous.  Thus,  the  objective  is  to  estimate  the  spectral  density  K(6)  = 
dF(9)/d9  and  we  restrict  our  attention  exclusively  to  estimators  defined  by 
rational  pr  functions  W(z). 

Let 

W{z)  -R 

=  W(zTTli" 

'  n 

denote  the  reflection  coefficient  of  the  load  normalized  to  R  ,  the  charac- 

n’ 

teristic  impedance  of  the  last  line.  Since  W(z)  is  pr,  p(z)  is  bounded-real 
(br);  i.  e. , 

|p(z)i<  1  ,  lz|  <  1  .  (82) 

Let 

2Q  (z) 
n 

denote  the  input  impedance  Z(z)  with  W{z)  =  R^,  Clearly,  this  corresponds 
to  p(z)  =  0  and  is  given  by  (74).  It  can  be  shown  (Appendix  A)  that 

Q  (z)  and  P^iz)  are  two  degree -n  polynomials  that  satisfy  the  even-part 
condition, 

It  is  also  a  consequence  of  the  same  analysis  that  the  even  part  of  the 
impedance  Z(z)  seen  looking  into  the  n+l-line  cascade  with  the  output  closed 
on  a  general  pr  termination  W(z)  is  given  by 

Z(z)  +  Z^(z)  l-p(z)(^(z) 

where 

D^(z)  =  PJ,z)-zp(z)P^(z) 

(17) 

A^(z)  =  A(l/z)  for  A(z)  real  and  rational. 


(8  5) 


(86) 
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1 


P^(z) 


zX^(z) 


(87) 


Consequently, 


Z(z)  generates  the  spectral  estimator 

l-|p(ej®)|^ 


K(0)  = 


Id  P 


(88) 


Observe  that  if  p(z)  sQ,  we  recover  Kj^j,(9),  Eq.  (73).  Clearly,  under  the 
assumption  that  F(6)  is  absolutely  continuous,  |  p(e'^®)  |  ^  1  and  W(z)  is  non- 
Foster.  In  general,  an  estimator  is  said  to  be  tuned  if  the  associated  re¬ 
flection  coefficient  p(z)  is  nonconstant. 

Let 


i.p(z)p^(z)  =  r(z)r,j.(2) 


(89) 


where  r(z)  is  a  real  rational  function  analytic  together  with  its  inverse  in 
|z|<  1.  (This  Wiener-Hopf  factor  is  uniquely  determined  by  the  normaliza¬ 
tion  condition  r(0)  >  0.  Furthermore,  r(0)  <  1  unless  p(z)  =  0.)  Hence, 
from  (85), 

K(e)=  |B(ej®)l^  (90) 


where 


B(.)  =  ^  .  (91) 

Since  D^(z)  is  free  of  zeros  in  |z|<  1,^^^^  B(z)  is  the  Wiener  factor  for  K(0)I 

Theorem  3.  Let  the  pr  function  W(z)  generate  the  spectral  estimator 
K(0),  given  the  error-free  covariance  data  C(k),  k=0-+n.  Let  {^j^}  denote 
the  sequence  of  radii  of  the  corresponding  circles  of  admissibility.  As  we 


' Pn(z)  is  the  polynomial  reciprocal  to  P^(z). 


(19) 


P  (z)  is  strict  Hurwitz  and 
n'  ' 


l-zp(z) 


n 


is  positive-real. 
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2 

know,  where  bQ  =  B(0).  In  addition,  from  (11)  and  (17), 

do  =  ^  /  inK(e)d0=  ZinbQ  .  (92) 

-rr 

Thus,  (91)  yields^^*^^ 

1)  =  r^0).^=r^0).r„(K^^j,)<r„(K^^)  (93) 

n-1 


with  equality  iff  K(9)  = 

2)  Entropy  (K)  =  entropy  (Kj^^,)  -in — ^ — <  Entropy  (K^^^)  (94) 


with  equality  iff  K( 0)  =  Kj^g(  0) .  Or,  as  an  alternative  to  (94)  we  may  also 
write. 

Entropy  (K)  =  entropy  (K^^,)  -  2^  /  in - T-g— d0  .  (95) 

.TT  ^  ■  I  n 


We  are  now  in  a  position  to  describe  an  entirely  new  type  of  spectral 
estimator  with  the  capability  of  controlling  both  resolution  and  robustness 
(as  measured  by  r^)  in  a  parametric  fashion.  To  this  end,  observe  that 
Eq.  (95)  permits  us  to  define 


• 


fn 


1 


p(e-'  ) 


j  03  I  2 


(96) 


to  be  the  entropy  echo-loss  density  at  angle 

A  flat-echo  estimator  (an  FEE)  is  one  for  which  the  entropy  echo-loss 
density  is  independent  of  9.  According  to  (96),  this  requires  that 

•  n 

lp(e'^”)|  =  constant  .  (97) 


But  this  is  possible  iff 

p(z)  =  ud(z)  (98) 

Bear  in  mind  that  Dj^(O)  =  Pj^(O). 

(21) 

'  '  The  phrase  "echo-loss"  is  taken  from  filter  theory  [7  ]. 
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where  [j  is  a  constant  constrained  to  lie  in  the  range 


0  <  n  <  1  ,  (99) 

and  d(z)  is  a  regular  all-pass;  i. e.,  d(z)  is  a  real-rational  function  analytic 
in  |z  I  <  1  which  satisfies  the  paraunitary  condition 

d(z)d^(z)  =  1  .  (100) 
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m.  THE  DESIGN  OF  TUNED  FEE'S 

Let  us  consider  the  properties  of  an  FEE  defined  by  the  reflection 
coefficient  p(z)  =  ud(z),  d(z)  a  regular  all-pass  and  0<  pi  <  1.  Then, 

r(z)  =  .  (101) 


^oo^*^FE^  -  (1  -  ^  ) 


■  00 


(102) 


and  Kpj.(9)  =  |B(e'^®)|^  where 


B(z)  = 


J  ^ 


Pj^(z)  -  Mzd(z)P^(z) 


(103) 


is  the  associate  Wiener-Hopf  factor.  Thus, 


1-U 


^FE<  =  SlE^  •  li.^ejVj®)U  (eJQ)!^ 


(104) 


where 


P„(2) 


n 


«„<">  =  in^r 

n 


(105) 


is  also  a  regular  all-pass  (because  Pj^(z)  Is  strict-Hurwitz).  Clearly,  Uj^(z) 
is  uniquely  determined  by  the  prescribed  data  C(k),  k  =  0-»n. 

We  can  now  list  some  of  the  anticipated  advantages  of  an  FEE. 

A^)  Numerical  robustness,  as  measured  by  r^,  is  assignable  in  ad¬ 
vance  by  fixing  the  constant  pi  and  is  unaffected  by  any  subsequent  tuning. 

A^)  On  lz|  =  1,  d(z)  and  Ujj(z)  are  both  of  unit  modulus  so  that 

d(e'^  ®)  =  explj  iji(9)]  (106) 

and 

U^(e->®)  =  exp  [j  9^(9)]  ,  (107) 
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(108) 


where  \|;{9)  and  are  real  functions  of  0.  Thus,  if 

cp(e)  =  t(0)  +  8^(0)+  e  , 

the  factor  that  multiplies  in  Eq.  (104)  reduces  to 

1-M^ 

Z”  • 

l-2ticoscp(0)  +  u 


Consequently, 


Kfe(9) 


l-|i  ^  ~ 


< 


i-m 

l-M 


(109) 


(110) 


and  any  such  ratio  can  be  realized  at  one  or  several  prescribed  values  of  0 
by  an  appropriate  choice  of  all -pass  d(z).  The  FEE  is  therefore  capable  of 
providing  selective  controlled  magnification  of  different  parts  of  the  spectrum. 
In  our  opinion,  if  used  properly,  this  estimator  could  prove  enormously 
valuable  in  situations  requiring  exceptional  resolution.  (The  design  of  d(z) 
is  one  of  phase  equalization  with  a  regular  all-pass  and  is  explained  in  detail 
in  Appendix  B . ) 

Aj)  The  poles  of  d(z)  are  zeros  of  B(z).  Hence,  the  FEE  not  only  in¬ 
cludes  the  MEE  as  a  limiting  case  (|j  =  0  =  maximum  robustness),  but  also 

yields  exact  model  matching  for  a  broader  class  of  random  processes  than 

(22) 

those  generated  by  pure  finite-order  autoregressive  schemes.'  ' 

A^)  Kf£(6)  >  0,  -TT<  0<  TT,  a  property  that  makes  good  physical  sense 
whenever  an  exact  spectral  null  is  highly  unlikely.  (Such  is  usually  the  case 
if  x^  is  the  sum  of  two  or  more  independent  random  processes.) 

Appendix  B  contains  design  formulas  for  a  singly  and  doubly-tuned  FEE. 
These  realize  maximum  gain  ( l+|j)/(  1 -ia),  at  one  and  two  prescribed  values  of 
0,  respectively. 

The  only  zero  of 

^  P  (z) 
n  ' 

is  z  =  00. 
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APPENDIX  A 

NETWORK  FOUNDATIONS 


Let  the  reactance  2-port  N^  shown  in  Fig.  1  be  described  by  its  2x2 
scattering  matrix  S(p)  normalized  to  on  the  right  and  R^  on  the  left.  Then 
[7  ]. 


S(p)  = 


g(z) 


h(z) 

g-npT 

Sii(p) 

Si2(p) 

-z"h^(z) 

_  ®12<P) 

"22<P) 

(Al) 


“  ZoT 

where  z  =  e  ^  and  h(z),  g(z)  are  two  real  polynomials  of  degree  n  in  z,  the 
latter  being  strict-Hurwitz.  Moreover, 

g(z)g,j.(2)  -h(z)h^(z)  =  1  (A2) 

Let  a  load  with  pr  driving -point  impedance  W(z)  close  terminals  ab. 
Then,  the  load  seen  from  terminals  ef  has  a  reflection  coefficient  normalized 
to  R^  given  by 

W{z)  -R^ 

s^  (z)  =  zp{z)  =  z  .  .  (A3) 


Hence  [12  ], 


s 


(A4) 


is  the  resultant  input  reflection  coefficient  normalized  to  Rq.  Making  use  of 
(Al)  and  (A2)  we  obtain, 

S(2)  =  MZ)  +.2P(Z)S(Z)  (AS) 

g(z)  +  zp(z)li(z) 


where  g  =  and  li  =  z^hjj, . 

Consequently,  the  input  impedance  is  given  by 
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where 


.p  l+sU)  . 

n  n 


_  g(2)  +  h(z) 

v'R^ 

p  (z)  =  . 


Note,  that  with  W(z)  =  P(z)  =0  and 


Z(z)  =  2j^£{z)  = 


2Q„(2) 


n 


Pn(-) 


Lastly,  by  a  direct  calculation. 


SO  that  (A6)  yields  (details  omitted), 

Z  +  Z^ 


1-PP* 


Thus, 


(P^-zpP  )(P^-zpP 
n  n  n  n 

Z(ej®)  +  Z^(e->®)  l-|p(ej®)|^ 


K(9)  - - j- 


in  which 


Dn^^)  =  ■  zP(z)P^(z) 


Naturally,  if  p(z)  =  0, 

«(«>  =  >S.E<«'  = 


(A6) 

(A7) 

(A8) 

{A9) 

(AlO) 

(All) 

(A12) 

(A13) 

(A  1 4) 


and  for  the  sake  of  completeness  we  shall  verify  the  determinantal  expression 
(74)  for  PJz). 
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From  (A9)  and  (AlO), 


^ME***'*‘^ME**“*  ,  '^n***^  _  1 

- 5 - ’TJiT^T^-pJ^rP-J^zV 


Let 


P  (z)  =  a.+a, z +'...  + a  ,  a->0  . 

n  '  0  1  n  0 


<A15) 


(A16) 


Then,  since 


Zj^j,(z)  =  C{0)+  ^2C(k)z‘"+0(z"+^)  , 

k=l 


(A15)  and  (A16)  yield 


00  00 

(C(0)+^  C(k)z‘‘"+  2  C(k)z‘")(a^+a^_jZ+ 
k=l  k=l 


+aoz’^)=: 


^  • 


(A17) 


(A18) 


It  follows  easily  from  (AlO)  that  P^(z)  ^  0,  lzl<  1,  so  that  l/Pj^(z) 
admits  a  power  series  expansion 


bQ+bjZ  +  b^z  +  . . . 


whose  radius  of  convergence  is  greater  than  unity.  Of  course, 


b.  = 


0  -  P^(0)  -  a^ 


=  —  >  0 


(A19) 


(A20) 


Consequently,  we  may  compare  the  coefficients  of  1,  z,  z^,  . . . ,  on 


both  sides  of  (Al8)  to  obtain 


T  a  =  b  e 
n-  o-n 


where  the  Toeplitz  matrix  is  given  by  (60), 


and 


2  =  ^V^n-l . ^0^' 

Sn"  (°’  0,  •  • ' ,  1) 


(A21) 


(A22) 


(A23) 
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is  composed  of  n  O's  and  a  single  1  in  the  last  row.  Thus, 


P^(z)  =  (z 


n-1 


l)a  =bQ(z 


n  n-1 


'  n 


(A24) 


C(0) 

C(l) 

.  .  .  C(n) 

C(1) 

C(0) 

.  .  .  C(n.l) 

C(n-l) 

C(n-2) 

.  .  .  C(l) 

n 

z 

n-1 

z 

...  1 

(A25) 


However,  if  we  now  use  Cramer's  rule  to  solve  Eq,  (A21)  for  aQ  we 


obtain 


(A26) 


and  (A25)  immediately  reduces  to  (74),  Q.E.D. 
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APPENDIX  B 


PHASE  INTERPOLATION  WITH  A  REGULAR  ALL-PASS 

The  phase  of  a  regular  all-pass  can  be  made  to  assume  any  prescribed 
set  of  real  values,  mod  2tt,  at  given  values  of  6=  9^,  02#  • . . ,  0^  .  However, 
the  degree  of  this  all-pass  can  far  exceed  the  integer  i  unless  the  phases 
and  the  angles  0^,  i=  l-»f  ,  are  interrelated  in  a  special  way.  In  this  appendix 
we  investigate  the  possibility  of  designing  a  degree -i  regxilar  all-pass  to 
interpolate  at  precisely  i  points.  This  of  course  represents  the  most  economic 
general  situation. 

Any  real  rational  regular  all-pass  d^  (z)  of  degree  i>0  has  the  structure 

2  (z)  z^E  ^(z) 

^  ^  •  E^TzJ  " 

where  €  =  ±  1  and 

E^(z)  =  l  +  fjZ  +  ...  +  f^  _jZ^ -^  +  1^2^  (BE) 

is  a  real  strict-Hurwitz  polynomial  with  lowest-order  coefficient  normalized 
to  unity  and  f^  ^  0. 

Since  d.  (1)  =  e,  the  problem  of  most  economic  phase  interpolation  may 

^  i  6  1 -s 

be  posed  as  follows:  Let  d^  (e"'  )  =  exp[j4>|  (S)  ].  Then,  1(^(0)=— ^tt  and 


CD 

1'#  (6)  =  (0)  +  arg  - ^ 

(B3) 

Given  the  two  sequences 

0=  0n<  9,  <  . .  .<  9,  <  TT 

0  1  1 

(B4) 

and 

o 

M 

P 

O 

P 

• 

• 

• 

P 

(B5) 

where  aj^=  (0j^)  -  (0),  find  E^  (z)  such  that 
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£,(e 

- 

E,(e 

The  original  recursive  p>plane  solution  to  this  problem  is  due  to 
Rhodes  [7,  13]  and  will  repay  careful  study.  Our  own  z -plane  development 
of  his  approach  places  less  emphasis  on  synthesis -theoretic  ideas  and  has 
a  somewhat  different  motivation. 

For  every  r  =  0-»i  ,  let 


,  k=0-»f  .  (B6) 


E,(g)  =1  +  5; 


k=l 


(B7) 


denote  the  real  strict  Hurwitz  polynomial  of  degree  r  which  interpolates 
correctly  at  0  =  6^,  0j^,  0^,  Then,  according  to  (B6)  the  difference 


^r+1  -  C  rd  d  1 

^_.^_c(dr^^-dr) 


(B8) 


±j0. 

must  vanish  at  z=l  and  z=e  ,  i=  I-»r.  Thus,  invoking  reality. 


£j.^l(z)  £y(z)  ri^(l-z)  .^^n(z^-2z  cos  0^+1) 


(B9) 


n  a  constant, 
r 


Specifically,  by  setting  z=0  it  is  immediately  verified  that 


r 


_  f(r+l) 
'  r+1 


(BIO) 


Clearing  fractions, 

^  (z^-2z  cos  ©i +1)  (Bll) 

i=l 


so  that 


rtZ^r+l  ■*'r+2  r+l  ’’rtl  ,  2  ,  „  ,, 

— - — - - - ^(z  -2zco.9_^^jtl).A^^,U).(B12) 

r+1  r  r+1  r 
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A  simple  rearrangement  now  yields 

^r+l^^r+Z'*' ^r+l®r^  ""  ^r+l^^r+Z"'  ^r+l^r^  • 

Since  E^^^(z)  is  strict  Hurwitz,  it  is  relatively  prime  to  £^^^{z)  and  must 
therefore  divide  +  A^^^(z)E^{z)  without  remainder.  Clearly,  the 

quotient  q(z)  is  at  most  of  degree  one  and  self-inversive;  i.  e, ,  q(z)  =  q(z). 
Consequently, 

q(z)  =  .  (B14) 


5  .  ,  a  real  constant,  and 
’  r+1  ’ 


’^r+l  2 


Obviously, 


ti  _  (r) 

?  -  1  X  r+1  _  r+z  r 

“r+1  "  "nr  ~  £(r+l)  .  f(r) 


(B16) 


r+1  r 


and  (B15)  takes  the  final  convenient  form. 


-  ?  x+l^  +  ,  {B17) 


Eq(z)=  1,  Ej(z)=  1  + 


sin^(9^-ai) 

sin^2Te]+^ 


r  =  f  -2 


Clearly, 


3in  jOj-Oj) 
sin  j(0j+a^) 


and  E^(z)  is  strict  Hurwitz  iff 


(B18) 


(B19) 


(B20) 
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This  granted,  we  now  prove  that  the  remaining  polynomials  are  also  strict. 
Hurwitz  provided 

0<  1  ,  r  =  l->i  -1  .  (B21) 

The  proof  is  inductive  in  character  and  divides  naturally  into  several  parts. 
1)  From  {B17),  by  direct  calculation  and  iteration, 

£  -E  =(?  ^,-l)(2^-2zcose  ^,  +  l)(£  -E  ^,£  ) 

r+2  r+1  r+2  r+1  '  r+1  ''  r+1  ''  r+1  r  r+1  r' 

(B22) 

2 

=  'nr^j(l-z)*  n  (z  -2z  cos  6^+1)  (B23) 

1=1 


where 

-1)  •  n  (§-1)  ^  0  .  (B24) 

Thus,  Ej^(z)  has  exact  degree  k,  k=  0-»f  , 

2)  Let  us  assume  that  all  Ej^(2),  k=  l-»r+l,  are  strict-Hurwitz  and  let 
us  write  E^^2(*)  =  ^r+2^^’  ^r+1^  stress  its  dependence  on  the  choice  of 
parameter 

Clearly,  as  is  seen  from  (B17),  E^^^l^rO)  has  two  distinct  zeros 

^1,2  ■ 

on  Iz]  =  1,  while  the  remaining  r  lie  in  |z|  >  1.  As  is  increased,  z^  and 

z^  either  slide  together  on  |z|  =  1,  or  move  as  a  pair  into  |z|  <  1,  or  penetrate 

jointly  the  region  [zj  >  1.  In  the  first  case,  there  exists  0<  <  1. 

such  that  E  .  t(z,  §  .  has  a  zero  z  =  e^®  of  unit  modulus. 
r+2  r+1 

Since  1)  has  the  zero  z=-l  and  r+1  zeros  in  [z]  >  1,  it  is  obvious 

that  in  case  two,  at  least  one  root  locus  must  leave  the  interior  of  the  unit 
circle.  Hence,  again  there  exists  a  ^r+1  ^  ^r+2^^’ ^ r+ 1^ 

has  a  zero  z  =  of  unit  modulus. 

Consider  case  three  and  suppose  that  for  some  0<  ?  <  1,  %) 

is  not  strict-Hurwitz.  By  continuity,  this  can  only  be  realized  if  some  root 
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locus  located  in  [zj  >  1  crosses  the  boundary  jz]  =  1  for  some 

0<  summary,  the  assumption  that  ?)  is  not  strict- 

Hurwitz  for  some  §  in  0  <  §  <  1,  necessarily  implies  the  existence  of  a 

i  6 

? 0<  <  1.  for  y/hich  has  a  zero  z  =  e"*  of  unit  modulus. 

This  we  now  show  is  impossible. 

3)  If  z  =  e'^®  is  a  zero  of  it  is  also  a  zero  of  £^^2^2, 

and  according  to  (BE3)  and  (B24),  either  6  =  0  or  it  coincides  with  some  6^, 
i  =  l-+r+l. 

If  6=0,  the  zero  is  at  z=l  and  (B17)  yields 

However,  E^(0)  =  E^^^(O)  =  1  so  that  E^(l)  >  0  and  E^^j^(l)  >  0.  (If,  say, 

E^(l)<  0,  then  E^(z)  must  vanish  for  some  zin0<z<l,  a  contradiction.) 
Consequently,  the  right-hand  side  of  (B26)  is  a  sum  of  positive  terms  and 
we  have  an  inconsistency. 

If  6=  6^^j,  has  the  zero 

z  =  exp[j6^_j_^] 


and  (B  17)  gives 


(BZ7) 


But  cos(6^^j/2)  /O  because  0<  and  (B27)  is  therefore  invalid. 

Lastly,  suppose  that  for  some  i,  i  =  1-^r,  6=6^.  Then,  from  (B17), 

0=§r+i(l+e  ^)  +  2(l-?^^^)(cos  e.-cos  e^^j)E^{e  ')  .  (B28) 


By  construction  and  the  fact  that  i£  r, 


j9-  j9. 

e„,(.  ■>  £,(e  •) 

"* - Ttr  =  - 3T" 


(B29) 


'The  coefficient  of  the  highest  power  of  z  in  Er+2(^»  ?r+l)  is  unequal 


to  zero  for  ail  in  0<?j.+  l<f«  This  means  that  no  root  locus  can  pass 

from  the  exterior  to  the  interior  of  the  unit  circle,  or  vice-versa,  by  escaping 
through  the  point  at  infinity. 
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and  since  for  any  k 

argfj^(e-^®)  =  k0-argEj^(e'^®)  ,  (B30) 

we  deduce  from  (B29)  that 

j0.  0.  j0. 

argE^^j(e  ")  = -^  +  arg  E^(e  ')  .  (B31) 

With  the  aid  of  (B31),  Eq.  (B28)  is  now  reduced  to  the  purely  real  form. 


0. 

j®i 

j®i 

0=  ^  ,  cos  • 

r+1  2 

+  (l-?r+l)(‘=°®  9i-cos  0^^^)  . 

E^(e  ") 

(B32) 

But  a  moments  thought  reveals  that  we  are  once  again  faced  with  a 

contradiction  because  0<  0.  <  0  ,  ,  <  TT  implies  that , 

1  r+ 1 

cos  9^-  cos  9^^^  >  0  .  {B33) 

Thus,  it  has  finally  been  established  that  every  polynomial  E^(z),  r  =  0-*t  ,  is 
strict-Hurwitz,  Q.E.D, 

The  restrictions 


0<c  <1,  r=l-»f-l  , 
r  ’  ’ 


(B34) 


impose  constraints  on  the  admissible  interpolatory  pairs  (0^,(1^),  i  =  0-+i  , 
These  constraints  reflect  the  well-known  fact  that  the  frequency  variation 
of  a  regular  all-pass  is  not  entirely  arbitrary. 


E^(z) 


For  example,  as  a  consequence  of  the  strict-Hurwitz  character  of 
we  can  easily  prove  that 


d 


) 

) 


dt^(0) 


>  0 


(B35) 


and 


=  kTT 


(B36) 


Thus,  the  a's  must  at  least  form  a  monotone -increasing  sequence: 
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(B37) 


0  =  Oq  <  <  . . .  <  <  i  n  , 

Unfortunately,  {B37)  does  not  always  guarantee  {B34)  and  the  general 
situation  is  rather  complicated.  Nevertheless,  because  of  their  great 
practical  importance,  we  shall  carry  out  the  designs  for  a  single  and  doubly- 
tuned  FEE  in  detail.  Recall  that  according  to  Eqs.  (104),  (108)  and  (109), 


‘^FE<8I  =  'Sv1eI9I 


l-2licos  cp(0)  +  u 


(B38) 


where 


9(6)  =  H9)  +  e„(e)  +  6  , 


(B39) 


i  6  i  6 

and  '([(0),  phase  angles  of  d(e‘*  )  and  U^(e‘^  ),  respectively. 

Design  of  a  Singly-Tuned  FEE.  Suppose  that  maximum  gain 


(B40) 


is  desired  at  some  single  suspected  peak  observed  in  say,  at  location 

6=0^.  Clearly,  we  must  have  9(6^)  =  0,  mod  Ztt,  and  this  we  attempt  to 
achieve  with  a  regular  all-pass  d^(z)  of  degree  one. 

Let 


(B41) 


so  that  0<  Yj  <  2tt  and  let' 


e  =  sign(Y^  -  tt)  .  (B42) 

If  e  =  0  we  pick  d^(z)  =  -1.  Otherwise,  if  e  =  +  l  we  choose  arg  11^(91)  =  2tt-Y^ 
and  if  e=-l  we  set  ilJj(6j)=TT  +  aj=2TT-Yp  Thus,  from  (B 18), 


sign  X  = 


1, 

x>0 

0. 

x=  0 

-1, 

X 

A 

0 
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sin  y  (6.  +  Y  +  Ti ) 

E,(z)  =  1  +  - f - - - —7^ - -  (B43) 

sin  tt) 


and 


(B44) 


It  is  easily  confirmed  that 


sini(9^+Y^+-4^TT) 
sin  ■  Y2+  tt) 


(B45) 


has  magnitude  less  than  one,  as  required. 

Design  of  a  Doubly-Tuned  FEE.  Let  us  now  attempt  to  achieve  maxi¬ 
mum  gain  at  two  locations  6^  and  0<  9^  <  9^  <  TT,  with  a  regular  all-pass 
d2{z)  of  degree  two.  Let 

Yi  =  9n( ®i)  ®i  »  mod  2tt  ,  (B46) 

i=  1,2  . 

Since  ®  1  ~  ®2  since  the  phase  angle 

regular  all-pass  ?,^{z)/E.^(z)  must  be  monotone -increasing, 
it  is  necessary  that  a2<  0i2* 

This  latter  observation  leads  to  the  following  assignment, 

1)  IfYj>Y2:  e  =  1,  a^  =  2TT-Y^,  i=l,2. 

2)  If  Yj<  Y2  ,  Yj<  TT,  Y2  >  tt;  €=  -1,  a^  =  TT-Yj  ,  a2=  3TT-Y2. 

3)  Any  other  case  is  impossible  with  a  degree-two  regular  all-pass. 

In  cases  1)  and  2),  E^(z)  is  given  by  Eq.  (B43). 

To  calculate  ^^(z)  we  can  use  (B17)  with  r  =  0: 

E2(z)  =  ?  ^(l+z)Ei(z) +(l-?j)(z^- 2z  cos  9i  +  l)  .  (B47) 

— m - 

'The  maximum  phase  increment  attainable  with  d2(z)  over  0<  9<  TT  equals  2tt 
and  occurs  at  9=Tr,  ~ 
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Of  course,  is  unknown  and  must  be  determined  from  the  interpolatory 
condition,  -o 

£  (e 

arg - jg -  =0.^  .  (B48) 

Ejle  ^ 

Clearly, 

j®2  j®2 

E2(e  n=?j(l+e  ^)E^{e  ‘')  +  2(l-§^)e  ^(cos  G^-cos  0^) 

and  we  find  that 


y  +  X 
~  y  +  X 


where 


y=E  (e  ‘•)e 


j0  -j9,/2 


X  = 


1-5^  (cos  G^-cos  Gj^) 


(B49) 


(B50) 


COS 


Finally,  solving  for  %  ^  with  the  help  of  {B43)  we  get  (algebra  omitted), 


02 

T 


_  J  cos 

^1  "  ^  cos  9^-cos  9^ 


1-e. 


sin ^  i, 


(B51) 


Let  us  check  to  see  if  0<  1  when  e  =  l.  For  this  it  suffices  that 


sin|-(92+Y2)  sinj(0j  +  Yj) 


sinl(02-Y2)  sini(0j-Yj)^ 


$  0  , 


(B52) 


depending  on  whether 


h-hi  °  • 


(B53) 


respectively.  These  inequalities  are  not  automatic  and  we  shall  attempt 
to  clarify  this  matter  in  Part  II. 
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PART  II 

A  DIGITAL  IMPLEMENTATION 
by 

OM  Gupta 

Digicomp  Corporation  Ithica  NY 


Chapter  3 


1.  Introduction 


The  solution  to  the  problem  of  the  spectral  estimation  as  proposed 
by  Professor  Youla  of  BrocAlyn  Polytechnic  offered  many  attractive 
features.  One  of  them  is  the  ability  to  tune  the  spectrum  at  desired 
frequencies.  The  technique,  called  Flat  Echo  Estimator  (FEE)  is  described 
in  a  pc^)er  by  Professor  Youla. ^  A  ocnputer  program,  described  in  the 
follomjig,  was  developed  to  inplement  this  technique  in  FORIEAN  IV  on 
MULTICS.  The  brief  description  of  the  cdgorithm  is  followed  by  the  block 
diagram,  description  of  various  siisroutines  and  a  listing  of  the  FORTRAN 
program  itself.  A  user  docixnentation  to  run  this  program  is  also  provided 
in  this  report. 

2.  Algorithm 


Youla  has  shown  that  the  maxinun  entrc^  spectral  estimates  I^(d) , 
for  the  finite  sequence  of  error-free  covariance  sanples  C(l),  C(2),..., 
C(n+1)  given  by 


(1) 


can  be  viewed  as  the  real  part  of  the  driving-point  impodenoe  of  a  cascade 
of  (n+1)  ccnmensurate  lines  terminated  cxi  a  load  W(z)  (“R^)  •  resistance 
is  the  same  as  the  characteristic  inpedanoe  of  the  last  line  in  the 
cascade.  He  showed  that  any  other  positive-real  load  termination  will  also 
provide  an  admissible  spectral  estimate  for  the  same  covariance  sanple 
sequence.  If  the  load  network  is  an  all-pa;v  nei^wJ;  tvith  reflection 
factor 

p(e^®)  =  Md(ei®)=  yexp(jt|((e))  ,  ,  (2) 

tlic  spectral  estimates  provided  will  be  Kf^(8)  arid  could  be  derived  to 
be 


Kfe<®>  = 


1-u 


2 


l-2u  cos  ti;(0)  +iJ 


2 


(3) 


Vhere 


C(0) 

The  phase  factor  6j,(0) 
ejqj  (j  SpO)) 


=  i|((0)  +  6r){0)  +  0 

is  given  by  the  expression 

=  iLi£L_- , 

Pn(2)  Pf>(z) 


(4) 


IM 


Hie  polyncmial  Pjj(z)  here  is  the  same  as  that  in  the  equation  (1) 
and  is  given  by 


1  C(1) 

Ci2) 

...  C(n+1) 

1 

1 

1  1  C{2) 

C(l) 

•  •  •  C  (n) 

P^(z)  =  -  1 

•^^n^n-l  !  •  •  • 

1  C(n) 

1 

C{n-1) 

...  C(2) 

1 

...  1 

(5)  : 

j 

Hie  determinant  is  given  by 

1 

r  C{1) 

1 

C(2) 

...  C(k+1)‘ 

- 

1  C(2) 

C(l) 

...  C(k) 

A.  = 

k 

1 

(6) 

L  C(k+1) 

C(k) 

...  C(l)  1 

Hie  factor  G(0) 

G(0)  =  - 

1-21100S  i(»(0)  + 

m2 

(7) 

in  the  eiqiression  for  Kpg(e)  is  similcu:  to  a  tuned  anplifier  qsdn  factor. 
Hie  maxinun  antropy  thus  can  be  tuned  to  provide  for  peaks  at  any 
desired  frequency  by  a  proper  choice  of  the  all-pass  load  network's  phase 
dicuracteristics  >1/(9),  Hie  choice  of  (0Mj<1)  provides  for  adjusting  the 
amplitude  of  the  peak.  It  should  by  noted  that  reduces  to 
Professor  Youla  has  demonstrated  a  link  between  the  vcuriable  u  and  the 
nunerical  robustness  of  the  spectral  estimator.  Hie  condition  y=0 
reduces  Kp^  to  and  represents  the  maximun  numerically  robust  choice. 

3.  Conputer  Program 

Hie  flow  chart  to  oonpute  Kpg(e)  for  a  given  sequence  C(l),  C(2),..., 
C(n+1)  is  showi  in  Figure  1.  Hie  first  step  to  obtain  lW(e)  is  to  oonpute 
K|^(e) ,  the  spectral  estimates  by  the  maxinun  entropy  metraxi^  then  design 
an  all-pciss  network  and  finally  oonpute  G(e) ,  Hie  multiplication  of 
G(e)  and  ^^(6)  will  provide  the  value  of  Kpp.(e).  It  should  be  noted  that 
we  don't  have  to  find  the  actual  element  vcilues  of  the  2d.l-pass  networic. 


FIG.  1 


:  FLOW  CHART  TO  ESTIMATE  SPECTRUM 
BY  FLAT  ECHO  ESTIMATOR  METHOD 
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Ratlter,  the  oatputation  of  an  appropriate  all-pass  bounded-recd.  ftxxrticn 
d(z)  (see  equation  2)  will  be  enough. 


3.1 


Ccmputation  of 


“Ihe  K^t;(e),  the  spectral  estimte  by  maxinun  entropy  method  is 
ocrputed  cis  follows.  First  a  matrix  t^^  is  formed  from  the  oovarianoe 
ooefficients  C{1),  C(2) ,  ...  C(npts)  as 


C(l) 

C(2) 

...  C(npts) 

C(2) 

C(l) 

...  C  (npts-1) 

(8) 


C(npts-l)  ...  ...  C(2) 


Ihen  the  ooefficients  of  a  polyncznial  C^(z)  (which  is  the  sane  as  P|^(z) 
except  for  the  constant  multiplying  feKTtor  (  =  1//  A^Ajj_j^)in  ecp\,  (5)  ) 
are  obtained  as  the  determinant  of  a  matrix  t^  formed  by  deleting  an 
appropriate  oolutn  i  of  matrix  t^^.  ^ 


Ihe  coefficient  of  for  example^  is  found  ais  the  determinant  of  a 
matrix  tj^^  obtained  by  deleting  the  first  oolutn  of  t^.  The  reascn  for 
doing  so  can  be  seen  by  noting  that  Quiz)  is  the  deteminant  of  the 
matrix  formed  by  t^^  appended  by  a  vector  as  its  last 

row,  i.e.,  _ 


npts 


•^AnAn-1  Pn<2)  =  Qn^z)  =  L  ...  1 


npts 

=  r 

m=l 


(-1)™ 


(9) 


The  ocrputation  of  det  performed  by  calling  the  subroutine 

pme  oof. 

A  constant  multiplying  factor  ^^It  is  obtained  as  the  product  of 
det  n  (=Aj^)  and  .  It  can  be  seen  that  A^^j^  is  same  as  the  constant 

term  in  Q^(z) .  The  determinant  Aj^  is  obtciined  in  the  do  loop  to 
CGfipute  ooefficients  of  0^(2). 


Finally,  the  coefficients  "ooef”  of  polyncndal  are  obtained  from 

the  cx3efficients  of  (^(z)  by  dividing  each  one  of  them  by  1//Qnult, 

Hie  spectrum  K^(e)  may  be  obtaiined  at  this  point  as  - rr — ^ 

|Pn(e^  )\ 

if  so  desired.  The  flow  chart  to  ocrpute  I^^(0)  is  shown  in  Figure  2. 

3.2  Ccnputation  of  FTE  Spectral  Estimates 

Our  goal  is  to  obtain  Kpg(e)  •  Therefore,  instead  of  siitply  obtaining 
K^^(e)  at  the  above  step,  a  subroutine  spectredL  estimator  is  called. 

Hiis  subroutine  first  queries  if  only  is  desired  or  Kp^  is  desired. 

For  I^,  it  asks  further  questions  as  to  singly  tuned  or  doitoly  tuned, 
the  value  of  these  frequencies  (in  radians)  and  the  maximun  robustness 
factor  AMU.  It  then  calls  an  appropriate  subroutine  design-all-pass 
or  design-all-pass-2  to  determine  a  first  order  or  second  order  cdl-pass 
reflection  factor.  The  subroutine  "response"  is  used  to  otxnpute  the 
spectrvm  nres  points  for  6  between  0.0  and  "thmax".  Hxai  return  fmn 
spectral-estimtor,  tiie  results  can  be  written  in  an  output  file  to  provide 
for  recidouts  and/or  plotted  on  Tek  screen  by  calling  plot  routines 
appropriately.  Hie  flow  chart  of  the  ocmputer  programs  is  shown  in 
Figure  3. 

4.  Subroutine  Description 

The  ocmputer  program  has  been  broken  down  into  very  many  smaller 
Droqrams  called  siiiroutines.  A  list  of  subroutines  used  herein  is  shown 
with  their  interconnection  in  Figure  4.  A  brief  description  of  various 
subroutines  is  provided  in  the  following: 

4.1  Subroutine  "Spectral-Estimator" 

4.1.1  Purpose:  The  purpose  of  this  subroutine  is  to  query  the  user 
if  m:  or  MEM  spectral  estimates  are  required.  In  case  FEE  is  desired, 
it  queries  the  user  about  the  robustr«ss  factor  AMU  and  the  tuned 
frequencies.  It  then  designs  an  appropriate  all-pass  reflection  coefficient 
(first  order  for  singly  tuned  by  calling  "desiqn-all-pass"  or  second  order 
for  doubly  tuned  by  calling  "design-all-pass-2".  The  queries  for  FEE 
parameters  and  the  design  of  all-pass  network  is  skipped  if  only  ^CM  is 
desired.  It  then  calls  "response"  to  ocTT«te  the  spectral  estimates 

using  equation  (3) .  The  flow  diagram  of  this  sul)routine  is  shewn  in 
Figure  5. 

4.1.2  List  of  Variables:  Following  is  the  list  of  variables  used 
in  "Spectral  Estimator" 

nfoe  -  integer  variable  =  0  for  MEM 

=  1  for  singly  tuned  FEE 
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COMPUTE 

SPECTRAL 

ESTIMATES 


PLOT  THE 
SPECTRUM  ANO 
WRITE  THE 
VALUES  IN 
OUTPUT  FILE 


Q  STOP  ) 


FIG.  3:  OVERALL  PROGRAM-FLOW  DIAGRAM 


FIG.  4:  INTERCONNECTION  OF  VARIOUS  SUBROUTINES 
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FEE? 


NO 


FIG.  5:  FLOW  CHART  -  SUBROUTINE  SPECTRAL  ESTIMATOR 
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=  2  for  doi±>ly  tuned  FEE 
amu  >  the  nimerical  robustness  factor 

peak  -  a  vector  of  frequencies  (in  radians)  at  which  tuning  of  the  spectnm 
is  desired.  'Ihe  frequencies  axe  cun:ai:>ged  in  the  increasing  order. 
Ihe  smallest  desired  frequency  must  be  peak  (1) . 
nerror  -  a  flag  to  indicate  that  all-pass  network  functicxi  cannot  be 
obtained 

ipts-d  -  the  order  of  the  polyngmieil  in  the  all-pass  reflection  factor 

E,(2) 

d(z)  =  - 

E,(z) 

npts-d  =  2  for  singly  tuned  case 
=  3  for  doubly  tuned  Ccise 
cf-d  -  the  coefficients  of  the  polynomial 

4.2  Subroutine  Etesign  all  pass 

4.2.1  Purpose:  This  subroutine  is  to  design  a  first  order  all¬ 
pass  network  fimction  d(z)  vrfvosc  phase  is  prescribed  at  one  frequency 
(peak  (1) ) .  This  oomputer  program  is  based  on  the  algorithm  provided  by 
Youla  in  Referenoe  1.  First  the  phase  GRFtIA  1  =  (6(ei)+  6j^)mod  2if 
at  0^  =  peak  (1)  is  computed.  The  all-pass  pluee  function  is  to  be 
chosen  in  such  a  way  that  i|)(6i)  +  GRNMh  1  ==  [0]  mod2iT  where  i|»(6j^)  is 
tie  phase  of  the  all-pass  d(z)  at  peak  (1) . 


4 . 3  Subroutine  amp-arq 

4.3.1  Purpose:  Given  a  polynomial  R(z) 

R(z)  =  aQ  ♦  aj^(z)  +  a.2^^  +...+ 

This  subroutine  ocmputes  the  anpTitude  and  phase  or  R(a)  for  any  givei  "th" 
z  =  e^^^^ 

The  subroutine  polyval  is  used  to  ocrpute  the  complex  value  of  polynomial 
R(z)  and  argnnt  is  used  to  determine  the  phase  of  R(z)  between  0  and  2it. 


It  is  frequently  desired  in  this  program  to  eveduate  the  phase 
all-pass  function. 


d(z)  =  0^^ - 


Ej(1/2) 


of  the 


E,(z) 


Efc(z) 


arg  fd(z)]  =  y  +  (n-l)e  ^  arq  (E  (z*)l  -  arg  [E  (z)] 

=  V  +  (n-l)e  -  2a (6) 
and  Y  is  sane  constant  phase  (0  or  2ti) 
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4.4  Subroutine  Polyval 

Hiis  subroutine  oanxites  the  value  (ooiplex)  of  a  polynotiial. 

R(z)  =  aj^  ♦  a2Z  +  a^z^  +...+  a^z*^"^ 
for  oonplex  z.  Ihe  coefficients  a^  are  real. 

4.5  Subroutine  Argrinent 

This  subroutine  calculates  tfie  arqument  of  a  ccxnplex  nuttfcer.  Hie 
result  is  between  0  and  2ir  depending  iqxin  the  signs  of  the  real  and 
imagina^  parts. 

4.6  Subroutine  "response" 


This  subroutine  Ccilculatcs  the  spectral  estimates  for  or  Kpg 
using  appropriate  expressions  at  nres  points  for  0  ranging  from  0 
bo  "thmax".  For  maximun  entropy#  Kjije  is  given  by 


fx(i)  =  Kj^{0i) 


1 


Fbr  is  multiplied  by  a  factor  Amult  given  by  equation  (7) . 

Most  or  the  variables  used  in  tliis  routine  aure  self-explanatory. 


4.7  Subroutine  design-all-pass. 2: 


4.7.1  Design  of  a  doubly  tuned  FEE 

The  ixirpose  of  this  routine  is  to  design  an  all-pass  network 
function  vhose  phase  has  been  prescribed  at  two  frequencies.  Except 
for  sane  minor  changes,  the  algorithm  used  is  the  same  as  provided 
by  Youla  in  Referenced) .  The  first  step  in  the  process  is  to  cor^iutt^ 
the  desired  phases  Ganma  1  aind  Ganra  2.  (eq.  B(46)  of  Reference  (1)  ) 
at  peaik  (1)  and  peak  (2) .  Then  the  phases  alpha  1  and  alpha  2  are 
OGiiputed  depending  upon  the  relative  values  of  Ganma  1  and  Garma  2. 


Vfe  know  that  (see  eq.  B48  of  Reference  {!)) 
152(2)  I 

arq  -  1  . 

E2<=')  I  ,  j0 


(4.7.1) 


z^((E2(i/z))  jze- 

-  =  e  -2arg(E2(z)) 

£2(2) 


=  e 


II-ll 


arg(K2(z)  ] 


j*!- 


=  02  -  0*2/2 


(4.2.2) 


E=e 


But  also  (see  eq.  (B-47)  of  Reference  (1) ) 


je, 

^2^®  )  =C^(l+e  )  *  t:^(e  n  +  (1-C^)(e  -2e  0080^^  +  1) 


^62/2  j02  j02 

=  2c^(cos02/2  •  e  )Ej^(e  )  +  (1-Cje  (2oos02  -  2oos0j^) 

j02  j02  16^2  OOS02  -  cosOj^ 

=  2^^  COS  2/2  e  (Ej^(e  )  +  -  .  - - - - 

COS02/2 


arg  E2(z) 


^®2 

=  arg(e  (y+x))  . 

362 

z=>e 

where  y  =  Ej^Ce  )e 

1-^1  00302-0080^^ 


00302/2 


From  eq.  (4,7.2) 

or 

i.e., 

but 


®2  ”  “2*^^  =  02  +  arg  (y+x) 

arg (y+x)  =  02/2  =  a 

Re(y+x)  =  R  cosA 
Iin(y+x)  =  R  3inA 

j02 

Ej^(z)  =  1  +  fj^z  =  1  +  fj^e 

"“^^2^2  “^6^2  _  j02 

y  =  E^(z)e  =  e  +  fj_e 


(4.7.3) 

X4.7.4) 


(4.7.5) 
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y=  COS02/2  -  j  sine2/2  cxjs02/2  +  sin02/2 


y=  (1  +  f2^)oos02/2  +  sine  2/2  (4.7.6) 

Therefore  frctn  eq.  (4.7.6) 

X  +  (l+fj^)cos02/2  =  RoosA . (4.7.7) 

(f  j^-1)  sine2/2  =  RsinA  .............  (4.7.8) 


Sincx2  is  known  (beinq  the  coefficient  of  z  in  Ej^(z)  of  sinqly  tuned 
case),  R  can  be  detennined  from  eq.  (4.7.8).  This  R  is  then  substituted 
in  eq.  (4.7.7)  to  obtain  x  . 

Once  X  is  known,  can  be  computed  as,  from  eq.  (4.7.5) 

X  00302/2 
oose^-oosSj^ 


1 


Then  the  coefficients  of  polynomial  E_(z)  are  determined  by  eq.  B(47), 
proof  in  Itef  (1),  as 

£3(2)  =  f,i^(l+z)  (l+fj^(z))  +  (1-C^)  (7.^-?r.  «-.Hej+:) 

=  5j^(l+(L+fj^)  z  +  fjZ^J  +  (1-Cj^)  (l-2z  oose^^  +  2^) 

=  1  +  ztC^d+f^)  -  2  cos0j^(l-Ci))  +  z^(eifi  +  (l-^i)) 
4.4.2  List  of  Variables: 

The  variables  used  in  this  subroutine  are  fairly  self-ex|5lanatory. 

Vector  cf-d  =  coefficients  values  of  polymial  £2(2) 
npts-d  =  nuptxjr  of  coefficients  of  £2^^) 

(=3  for  doubly  tuned  case) 
ar  =  A  of  eq.  (4.7.6) 
br  =  B  of  equ.  (4,7,9 


11-13 


5. 


User  Instructions 


Ihe  ocnputer  program  is  iRplemented  in  IV  and  the  listing 

is  shewn  in  ;^3pendix  A.  Hie  program  is  named  as  F1SE2  on  Multics. 

Before  running  the  program,  the  user  must  make  sure  that  the  gre^iiics 
mode  for  his  system  is  turned  on,  A  typical  user  session  is  described 
below.  Hie  system's  response  or  questions  are  marked  as  "S”  vAiereas  the 
user  response  is  marked  aa  U,  (see  figure  7) 

S:  Y 
U:  fee2 

S;  enter  two  digit  file  nixrber  (note  1) 
enter  number  of  parts  to  be  used 

S:  enter  0  for  mem,  1  for  singly  tuned  and  2  for  doubly  tuned  fee 
(note  2) 

U:  0,1,  or  2 

Hie  following  questions  are  not  asked  it  the  user  enters  0  in  the  above 
response. 

S;  enter  robustness  factor  -  amu 
U:  any  real  nvmber  u,  0<v<l  (note  3) 

S;  enter  desired  resonant  peaks  (in  radians) 

U;  any  real  nuntier  between  0  and  1  (note  4) 

S':  Hie  estimated  spectral  plots  are  plotted,  the  u  and  resonant  peak 
values  printed  and  then,  the  following  statement  typed:  If  you  wish 
to  make  a  hard  copy,  do  it  new. 

Hien  type  1,  if  you  wish  to  plot  the  phase  of  the  all  pass  network. 
Else,  type  in  0. 

U:  0  or  1 

S:  depending  upon  the  user's  response  to  the  above  question,  the  phase 
of  the  all  pass  FHE  network  load  is  plotted 
S:  Stop 
S;  y 

5.1  Notes: 

5.1.1  Note  1:  Hie  program  at  this  stage  is  asking  you  to  input 
two  digits  XX  (=25  in  test  case)  of  the  file  XX  on  the  terminal.  Hie  . 
program  assumes  that  the  following  input  data  is  pre-stored  in  file  XX 
(file  25, for  exanple) . 

Input  data  in  file  XX: 

(i)  nres  (nuriier  of  points  at  vt4iich  estimated  spectral  response 
is  to  bo  evaluated  in  format  14) 
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(ii)  thnax  (the  maxiimm  value  of  e  (in  radians)  up  to  \*^iich 

spectral  estimates  cire  to  be  oorputed  in  format  EIS^IO) 

(iii)  npts  (in  format  12,  the  number  of  error-free  covariance  sanple 

data  points) 

(iv)  Hie  oovcurianoe  values  at  npts  points  (in  the  format  E18.10) 

Bor  the  test  case,  the  data  stored  in  file  25  is  shown  in  fig. 6.  Hie 
MM  and  FEE  spectral  estimates  are  shown  in  figures  8  through  11. 

5.1.2  Note  2:  Hie  program  at  this  point  is  asking  if  the  maxinun 
entropy  spectral  estimates  or  FEE  estimates  are  desired.  If  the  user 
types  in  0  (for  MEM) ,  the  program  does  not  ask  any  further  questions, 
and  sinply  plots  the  MEM  spectrcd  estimates  on  the  terminal  screen.  In 
the  case  of  FEE  (user  response  of  1  for  singly  tuned  and  2  for  doubly 
tuned) ,  the  user  has  to  type  in  the  robustness  factor  w  and  the  desired 
resonant  peaks. 

5.1.3  Note  3;  Hie  robustness  factor  jj  should  he  a  real  nutber,  between 
0.0  and  1.0.  If  you  type  0.0,  the  vector  will  be  the  saome  as  MIM 
estimates.  Hie  program  would  bomb  if  u  =  1,0. 

5.1.4  Note  4;  Hie  values  of  resonant  peaks  in  radians  that  you  are  oslced 
to  input  (real  miibers)  must  be  1  for  singly  tuned  and  2  for  doubly 
tuned  case.  FXirthermore,  the  peaks  for  doubly  tuned  case  must  be  in 
increasing  order,  i.e.  the  first  input  value  must  be  less  than  the  second 
input  value. 

5.1.5  Note  5:  At  this  point,  the  system  has  oenputed  and  plotted 
the  desired  spectral  estimates.  It  has  cdso  printed  the  other  useful 
information,  such  as  w  and  peak  values  for  the  hard  copy  records.  Hie 
program  can  also  plot  the  phase  of  the  all  pass  function  d(Z) .  It  asks 
the  user  if  he  wishes  to  plot  this  phcise.  Since  the  program  shall 
blcink  out  the  screen  before  plotting  this  cdl  pass  phase  function,  the 
user  must  make  a  hard  copy  of  spiectral  estimates  before  respxmding  to 
this  question. 
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Fig.  6:  Data  for  Test  Exanple  (ocntained  in  file  25) 
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Setip  graphics-table  tek  4014 
segm^t  Setup-graphics  not  found, 
r  1217  0.055  1.750  50 

Setup  graphics-table  tek  4014 
r  1213  0.351  3.906  62 

fee2 

enter  two  digit  file  ntinber 
25 

enter  0  for  nen,  1  for  singly  and  2  for  doti^ly  tuned  fee 
2 

enter  robustness  factor  -:  anu 
0.5 

enter  desired  resonant  peak  (in  radians) 

0.6035 

enter  desired  resonant  peak  (in  radians) 

2.502 


Fig.  7:  A  typical  user  session 
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Part  I 


1.60p 


diinsnsion  c(20),  tn(20,20) ,  peak(20 
dimension  ooef  (20),  fx{1024) ,  thi(1024) 
character* 16  header  (3) 
ocnplex  z,  fz,  cexp,  cth 
oatmQn/respos/ph__edlpciss  (1024 ) 
external  plot^seti^j  (descriptors) 
external  plot_(descriptors) 
external  ioa_(descriptors) 
extenuil  read__list_(descriptors) 
extemcil  write_list_(descriptors) 
header  (1)  -  "Maximun  Entropy" 
header  (2)  -  "Singly  tuned  FEE" 
header  (3)  -  "Doubly  tuned  FEE" 
call  ioa_(enter  two  digit  file  nurtoer) 
call  read  list_(nfil) 
read(nfil,113)  nres 
read (nf il , 111) thmax 

113  format  (14) 

read (nfi 1,110)  npts 
do  10  i=l,  npts 

10  read (nfi 1,111)  c(i) 

110  format  (12) 

111  format  (E18. 10) 
nptsi=npts-l 

do  40  jr=l,nptsi 
do  20  i=l,jr 

20  tn(jr,i)  =  c(jrfl-i) 

kd  =  npts-jr+1 
do  30  i=2,kd 
jc=jr+i-l 

30  tn(jr,jc)=c(i) 

40  continue 

nfilot=18 
nsign=l 
deln=0.0 
do  60  i=l,npts 
call  pTEOof (nptsi,i,tn,cf) 
coef (npts-i+l)«nsign*cf 
dsin  =  (nsign*c(npts-i+l) *cf)+deln 
nsign  =  -nsiqn 

Appendix  -  Fortran  Program  Listing 
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(i^jpendix,  cxxitinued) 


60  continue 

write (nf ilot , 114 ) 

114  fonnat(Sx, "Spectral  Estimates  follow") 

anult  =  deln  *  ooef(l) 

anult  =  abs (anult) 
anult  =  1.0  /  anult 
anult  =  sqrt  (anult) 
do  75  i=l,  npts 

75  ooef(i)  =  anult  *  ooef(i) 

call  spectral_  estimator(thi,£x,coef,thmaX/anu,pea)c, 
npts ,nres ,n£es) 
do  80  i=l,nres 

write (nfilot, 115)  i,  thi(i),  fx(i),  ph_allpass (i) 

115  forTnat(I4,  5x,  E12.4,  4(2x,E14.5)) 

80  continue 

call  plot_Ssetup(header(nfee+l) ,"freq", "spectral 
value", 3, 10. 0,1,0) 
call  plot_(thi,fx, 1024,1,) 
if (nfee.Qq.O)go  to  95 

call  ioa__( "robustness  factor  and  peak  freq  are") 
call  write_list_(amu) 
do  90  i=l,nfee 

90  call  write_list_(peak(i) ) 

call  ioaC'Input  data  is  in  file  No.") 
call  write_list_(nfil) 

95  continue 

call  ioa_("If  you  wish  to  make  a  hard  copy,  do  it  now") 

call  ioa_("'Itien  type  i  if  you  wish  to  plot") 

call  ioa_("the  phase  of  all  pass  fee  network") 

call  ioa_("else  type  in  0") 

call  read_list_(ntenf3) 

if(nteiTf5  .ne.  1)  go  to  96 

call  plot_Ssetup( "phase  of  tuned  foe  network", "freq", 

"pliase", 1,10. 0,1,0  /c) 

call  plot_(thi,ph__allpass, 1024,1,"  ") 

96  continue 
close (nf il) 
close (nf ilot) 
stop 

end 
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(Appendix,  cxxitinued) 


subroutine  spectraI_estiinator(thi,£x,cf ,ttTiiiax,c3nu 
peak,npts,nres ,nfee) 

dirension  thi(1024)  ,fx(1024)  ,peak(20)  ,cf  (20)  ,cf-d(20) 

external  ioa- (descriptors) 

external  read-list- (descriptors) 

call  ioo- ("enter  0  for  mem,  1  for  sinqly  and  2  for 

doiisly  tuned  fee") 

call  read-list- (nfee) 

if (nfee.eq.O)  qo  to  20 

call  ioar  ("enter  rciwstness  factor  -:  ami") 
call  read-list- (amu) 
do  10  i=l,nfee 

call  ioa-("<aiter  desired  resonant  pea}c(in.itt>U<ms) "] 
call  read-list- (peaJc  (i) ) 

10  continue 

call  design-allpass  (npts,npts-d,cf,peaJc,cphase-d, 
cf-d) 

if (nfce.eq.l)  go  to  20 
if  (pea)c(2)  .gt.peaJc(li)  go  to  15 

call  ioa- ("error:  2nd  res  freq  less  than  1st  res  freq. 
Treated  as  singl/cy  tuned") 
go  to  20 

15  nerror  =  0 

ccill  desiqn-edlpass-2(npts,npt^  d,nerror,df,pea]c, 

<73hase-d,cf-d) 

if (nerror  .eq.  0)  go  to  20 

Ccill  ioa- ("error:  All  pass  can  not  be  designed.  Treated 
as  singly  tuned/c") 

20  continue 

call  response  (nres,nfee,npts,:^its-d,cf  ,thimax,cf-d, 

cphase-d ,  amu ,  &c ,  thi ) 

return 

aid 

subroutine  response(nres,nfee,npts,npts-d,cf,thmax, 

cf-d,ophaserd,ar!u,fx/c,thi) 

oarnnon/respos/ph-cillpass  ( 1024 ) 

nfee  =  0  if  mem  is  desired,  »  nirber  of  resonant 

pea}(S  for  fee) 

dimension  fx(1024),  thi (1024) ,c£(20) , cf-d (20) 

delth=nrcs 

delthothmax/delth 

anusq^anu^amu 

anult-nun  s  l.OefOO  -  cnusq 
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(A|)pctxlix,  continued) 


amlt-den  =  l.Oe+00  +ainLisq 
aniii2  =  anu+oiru 
do  80  i=l,nres 
th  =  i*delth 
thi(i)  =  th 
cphase  =  O.OetOO 

20  continue 

a  =  (peak(l)-alpiial)/2.0e+00 
b  =  (peak(l)+alphal)/2.0et00 
cf-d(2)  =  sin  (a) /sin  (b) 
cf-d(l)  =  l.efOO 
npts-d  =  2 
do  30  i»3,20 

30  cf-d(i)  »  O.efOO 

c^dl  anp-cirq  (npts-d, cf-d, peak (i)  ,cphcise-d,anT>,phase-d, 

all-pfiase-d) 

nfil=40 

101  fonnat{2x,  "Design",  4x5E13, 6) 

th  =0.eH)0 
close (nfil) 
return 
end 

function  twopi-inod(a) 
twopi  =  6,283184et00 
na  =  a/tMopi 
twopi-nod=-na*twt^i  +a 

if  (twopi-ftod. It, O.efOO)  twopi^noe^twopiftMopi-inod 

return 

end 

subroutine  design-eiLl-pass-2(npts,npts-d,nerror, 
cf ,peak,cphase-d,cf-d) 
diroensicn  cf  (20)  ,cf-d(20)  ,peak(20) 
c. .  .Ihis  subroutine  assunes  that  singly  tuned  beai  cedled 
cnoe 

pi  =  3.141592ef00 
twopi  =  6,283184e+00 
threepi  =  9,424696ef00 
cphase  =  O.efOO 

call  anp-arg(npts,cf,peak(l) , cphase, anpl,phasel, 
all-pha»3l) 

call  anp-arg(npts,cf  ,peak(2)  ,cphase,an^,phase2, 
all-phase2) 

gann^  =  all-phcisel  f  peak(l) 


A- 4 
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garina2«  all-phase2  +  peak  (2) 
garinal  =  twopi-modCganinal) 
ganina2  =tMDpi*4nod(ganina2) 
nfil=40 

call  aii|>-arg(npts-d,cf-d,peak(2)  ,cphase-d,app, 
ph2,aph2) 

call  anp-arg(rpts--d,cf-d,peak(l)  ,cphase-d,cil,phl,aphl) 
alphal  =  twcpi  -  ganinal 
alphci2  =  tMopi  -  ganina2 
eps  =  O.Oe  +  00 

if  (ganinal  .ge.  ganina2)  go  to  10 
if  (ganinal  .gt.  pi)  go  to  100 
if(ganina2  .It.  pi)  go  to  100 
eps  =  l.Oe+00 
alphal  =  pi  -  ganina  1 
£ilpha2  =  thre^i  -  ganina2 
10  oontiniie 

nfil=40 

ar  =  -0.5e+00  *  cilpha2 

r  =  (cf-d(2)-l.eH)0)*sin(peak(2)/2.0e+00)/sin(ar) 

X  =  (cf-ct(?)+1.0e+00)*cos(peak(2)/2.0e+00) 

X  =  -X  -  (r*oos(ar)) 

br  =  x*o.>8(peak{2)/2.0e+00)/(cos(peak(2))-oos(peak(l))) 

epsi  =  l.e+00/(l.e+00  +  br) 

epsmi  =  l.OefOO  -  epsi 

cf-d(3)  =  epsmi  +  (epsi*cf-d(2)  ) 

cf-d(2)  =  {epsi*(l.e+00+cf-d(2)))  -  2.0eM)0*c!Os(poak(l)) 
*^>snd 

cphase-d  =  eps  *  pi 
ipts-d  =  3 

call  aitp-arg(npts-d,cf-d,peak(2)  ,cpheBe-d,anp,ph2d, 
aph2d) 
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ABSTRACT 


Let  Xj  denote  a  discrete-time,  zero-mean,  second-order  stationary, 
full-rank  m-dimensional  random  vector  process.  Let  us  assume  that  the 
spectrum  of  x^  is  absolutely  continuous  and  let  K(6)  and  C(k)  denote  the 
associated  mxm  spectral  density  and  covariance  matrices,  respectively. 

Suppose  that  the  partial  information  C(0),  C(l),  •••,  C(n)  is  known  in 
advance  without  error.  The  first  major  contribution  of  this  research  report 
is  contained  in  theorem  2  which  supplies,  via  Eqs.  (173)  and  (174),  a  general 
explicit  parametric  formula  for  the  class  of  all  spectral  densities  K(0)  con¬ 
sistent  with  the  given  n+1  (matrix)  pieces  of  covariance  data. 

The  free  parameter  is  an  mxm  matrix  function  p(z)  which  is  analytic 
in  jz]  <  1  and  of  euclidean  norm  less  than  unity  almot  everywhere  on  [z]  =  1, 
Thus,  p(z)  is  essentially  an  arbitrary  passive  scattering  matrix| 

Most  of  the  numerical  calculations  can  be  carried  out  efficiently  and 
recursively.  Indeed,  our  algorithm  can  be  constructed  around  the  generalized 
Levinson  recursions  for  the  left- right  matrix  orthogonal  polynomials  generated 
by  the  sequence  C(0),  C(l),  ...,  C(n).  Consequently,  all  the  highly  developed 
and  sophisticated  software  that  has  come  into  existence  during  the  last  decade 
because  of  the  intense  effort  devoted  to  Burg's  maximum -entropy  estimator, 
can  be  used  intact. 

The  second  major  contribution  is  the  identification  of  the  FEE  or  Flat- 
Echo  Estimator.  The  FEE  is  based  squarely  on  frequency-domain  concepts 
and  represents  an  entirely  new  particular  solution  to  the  problem  of  spectrum 
recovery.  It  emerges  by  an  appropriate  restriction  of  the  functional  form  of 
the  parameter  p(z).  Some  of  its  principal  advantages  are  the  following. 

1)  "Numerical  robustness"  can  be  assigned  in  advance  by  the  adjustment 
of  m^  numerical  parameters. 

2)  It  is  possible  to  "tune"  the  estimator  to  produce  selective  spectral 
amplification  at  a  finite  number  of  desired  frequencies  without  impairing 
either  its  interpolatory  character  or  its  numerical  robustness.  Enhanced 
resolution  is  to  be  expected. 


1 

1 


3)  For  the  setting  of  the  numerical  parameters  which  yields  maximum 
robustness,  namely  p(z)  =  0^,  the  FEE  coincides  with  the  MEE. 

In  conclusion,  we  should  like  to  mention  that  our  development  of  the 
general  parametric  formula  for  K(3}  alluded  to  above,  constitutes  the  solu¬ 
tion  to  a  long  open  problem  which  some  researchers.  Burg  included,  have 
considered  to  be  prohibitively  difficult. 
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I.  INTRODUCTION 


In  a  previous  publication  [l  ]  the  author  developed  an  entirely  new 
solution  to  the  problem  of  spectrum  recovery  from  a  given  set  of  covariance 
samples  of  a  single  random  process. 

The  flat-echo  estimator  (the  FEE)  incorporates  the  following  important 
features : 

1)  Numerical  robustness  is  assignable  in  advance  by  the  adjustment 
of  a  scalar  parameter  and  can  be  traded  off  against  resolution. 

2)  It  is  possible  to  "tune”  the  estimator  to  provide  selective  amplifi¬ 
cation  at  any  set  of  prescribed  frequencies  without  impairing  either  its 
interpolatory  character  or  its  numerical  robustness. 

3)  For  the  parameter  setting  which  yields  maximum  robustness,  the 
FEE  coincides  with  the  MEE,  the  maximum -entropy  estimator. 

Our  principal  objective  in  this  continuing  study  is  to  generalize  the 
FEE  to  the  multichannel  case  and  we  therefore  assume  that  the  reader  is 
familiar  with  the  geometric  and  analytic  points  of  view  adopted  in  Ref,  I, 
Moreover,  since  this  is  a  research  report,  our  demands  on  his  knowledge 
will  vary  widely  but  we  have  nonetheless  attempted  to  be  tutorial  whenever 
possible. 

With  regard  to  the  matter  of  notation,  most  of  it  is  standard  and  self- 

_  /  *  _ / 

explanatory.  However,  let  us  point  out  that  A,  A  ,  A  (=  A  )  and  det  A  denote, 
in  the  same  order,  the  complex  conjugate,  transpose,  adjoint  and  determinant 
of  the  matrix  A. 

Column-vectors  are  written  a,  x,  etc.,  or  as 

X  =  (Xj,X2,  •  •  •  ,x^)'  (0.0) 

whenever  it  is  desirable  to  exhibit  the  components.  Further,  1  is  the  mxm 

m 

identity,  0  ,  the  mxk  zero  matrix  and  0  is  the  m -dimensional  column- 

^  m, k  -m 

vector. 

9j( 

As  usual,  if  A=  A  is  hermitean,  A>0(>0)  means  that  x#Ax>0(>0)  for 
every  vector  x  /  O.  We  then  say  that  A  is  nonnegative  or  positive -definite, 
as  the  case  may  be. 


In  this  paper,  we  distinguish  carefully  between  the  right-half  p-plane 
and  the  interior  of  the  unit  circle  since  some  of  our  arguments  are  couched 
in  both  analogue  and  digital  terms. 

Thus,  if  W(p,  z)  denotes  a  matrix  whose  entries  are  meromorphic 

:0 

functions  of  both  p=  O+jU)  and  z  =  r  e'^  ,  it  is  convenient,  for  typographical 
reasons,  to  set 


W^(p,z)  s  W*(-p,  1/z)  .  (0.2) 

Clearly,  for  p=juu  and  z  =  ,  tu  and  9  real, 

W^(jU),e-’®)  =  W*(j(JU.ej®)  (0.4) 

and  the  sub-star  coincides  with  the  top-star. 

Let  X  stand  for  either  p  or  z  and  let/? denote  either  the  open  right-half 
p-plane,  Rep>  0,  or  the  open  unit-circle,  jz]  <1. 

An  mxm  matrix  Z(x)  is  said  to  be  positive  if  1)  it  is  analytic  in  /?and 

2). 

>  0  ,  xe/e  .  (0.6) 

Similarly,  an  mxm  matrix  W(x)  is  said  to  be  bounded  if  1)  it  is  analytic  in/? 
and  2), 

1  -W(x)W*(x)>0  ,  xe/e  .  (0.8) 

m  — 

If  in  addition,  Z(x)  or  W(x)  is  real  for  all  real  xe/?,  it  is  said  to  be 
positive -real  or  bounded-real,  respectively. 

A  positive  matrix  Z(x)  which  is  a.e.  skew-hermitean  on  the  boundary 
of/?is  said  to  be  lossless  while  a  bounded  matrix  W(x)  which  is  a.e.  unitary 
on  the  boundary  of /?is  said  to  be  regular  paraconjugate  unitary.  If  Z(x)  is 
also  real,  it  is  Foster  and  if  W(x)  is  real  it  is  regular  paraunitary. 
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U.  FORMULATION 

Z  Let  denote  a  discrete-time  m -dimensional  random  vector  process 
where  "t"  can  traverse  all  positive  and  negative  integers.  Let  us  also  sup¬ 
pose  that  the  process  is  zero-mean  and  second-order  stationary.  Then, 

'  11^1=  0-x*  ,  (1) 

is  the  associated  mxm  covariance  function.  ^  Clearly,  for  all  integers  k, 

C(-k)  =  E(x^x*_j^)  =  =  C*{k)  .  (2) 

As  is  well-known  [2],  if  the  power  spectrum  F(6)  of  the  process  is 

absolutely  continuous,  there  exists  an  mxm  hermitean  nonnegative -definite 

2 

spectral  density  matrix  K(9),  whose  entries  are  absolutely  integrable  over 
-TT<  0<  TT  and  is  such  that, 

C{k)  =  ^  /^e’-'‘"®K(e)de  .  (3) 

-TT 

Hence, 

go 

K(e)-y  C(k)ej‘"®  (4) 

k=-OC 

and  the  mxm  covariance  samples  C(k),  lk|  =0-*oo,  emerge  as  the  Fourier 
coefficient  of  K(6). 

Moreover,  if  the  process  is  free  of  essentially  determined  components, 
K(0)  satisfies  the  Paley-Wiener  criterion, 

J  fn  det  K(  0  )d6  >  -  go  .  (5) 

-TT 

Very  often,  we  do  not  distinguish  between  the  functi on  C(k)  and  the 
particular  sample  C(k).  Hopefully,  when  taken  in  context  this  should  never 
cause  confusion. 

^For  short,  K(0)eLj. 
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Clearly,  such  a  process  automatically  obeys  the  full-rank  condition, 


det  K{  9)  ^  0  ,  a.  e.  , 


but  not  conversely. 

From  the  inequality  K(0)>  0  it  is  easily  shown  [6]  that 


(5a) 


'c(O) 

0(1) 

.  .  .  0(n) 

C(-1) 

0(0) 

.  .  .  C(n-l) 

_C(-n) 

C(-n+l) 

.  .  .  C(0) 

(5b) 


Furthermore,  if  the  process  is  full-rank  (and,  a  fortiori,  if  it  satisfies 
Paley-Wiener),  only  the  inequality  sign  prevails. 

Unless  stated  explicitly  otherwise,  we  shall  adhere  to  two  main  working 
assumptions . 

A^)  The  power  spectrxxm  F(0)  of  the  process  is  absolutely  continuous, 


i.  e.  , 

dF(0)  =  K(9)de 

(6) 

where  K(9)eLj  and 

K(9)  =  IC’-(6)  >  0  . 

(7) 

A^)  Either  the  process  is  full-rank  or  its  satisfies  the  Paley-Wiener 
condition  (5)  or  both.  The  precise  situation  can  always  be  ascertained  from 
the  context. 

Our  primary  task  is  to  develop  a  physically  insightful  solution  to  the 
problem  of  estimating  K(@)  from  a  knowledge  of  n+1  given  error-free  mxm 
covariance  samples  C(0),  C(l),  ...,  C(n), 
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ni.  THE  2m -PORT  INTERPOLATORY  CASCADE 


By  way  of  motivation,  let  us  recall  the  properties  of  a  cascade  N  of 
n  ideal  TEM  lines  which  have  the  same  1-way  delay  T>0  but  arbitrary  positive 
characteristic  impedances  Rq,  Rj^,  ....  R^^  ^  [l  ]. 

Let  Sj^(p)  denote  the  scattering  matrix  of  N  normalized  to  R^^  at  the  in¬ 
put  and  any  positive  port  number  at  the  output.  Then,  Sj^(p)  is  completely 

parametrized  by  the  specification  of  a  single  arbitrary  real  polynomial  h(z) 

3 

of  the  form 

h(z)  =  zh^+z^h^i- .  •  •  +  2*^^  .  (8) 


Explicitly, 


C^)  with  h(z)  so  prescribed,  choose  the  real  polynomial  g(z)  (of  degree 
<  n-1)  to  be  the  strict-Hurwitz  solution  of  the  equation 

g(z)g.j,(z)  =  H-h(z)lx^(z)  .  (9) 

This  Wiener -Hopf  factor  is  analytic  together  with  its  inverse  in  |z|  <1 
and  is  uniquely  determined  under  the  restriction  g(0)>0. 


C^)  Let 


Then, 


-2pT 

z  =  e  . 


h(z) 

g-npT 

e-npT 

-z'^h^(z) 

(9a) 


(10) 


It  is  readily  seen  that  Sj^(p)  is  a  real  meromorphic  matrix  function  of 
p  that  is  analytic  in  Rep>0  and  satisfies  the  identity 

Thus,  Sj^(p)  is  regular  paraunitary. 


The  property  h(0)  =  0  reflects  the  fact  that  the  input-side  normalization 
number  has  been  chosen  equal  to  Rq,  the  characteristic  impedance  of  the 
first  line. 
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The  immediate  goal  is  to  generalize  the  above  construction  to  the  2m- 
port  case,  m  >  1.  However,  due  to  the  lack  of  commutativity  of  matrix 
multiplication,  several  subtleties  arise  for  m  >  1  which  car.  only  be  dealt  with 
by  modifying  the  procedure  in  a  substantial  way. 

Lemm a  1 .  Let  H(z)  denote  an  arbitrary  (and  not  necessarily  real)  mxm 
polynomial  matrix  of  degree  <  n  of  the  form 

H(z)  =  zH,+z^H,+  -- •  + z”h  .  (12) 

'12  n 

There  exist  two  mxm  polynomial  matrices  Gj(z)  and  0^(2)  of  degrees  <  n-1 
such  that 

1)  det  G^(z)  •  detG2(z)  j^O,  iz|<l.  (13) 

2)  The  matrix 

L(z)  =  G"\z)H(z)G2(z)  (14) 

is  polynomial. 


is  meromorphic  in  p,  analytic  in  Rep>0  and  paraconjugate  unitary: 

Sn(P>Sn<.Ip>  =  hm  ■ 


Proof.  Let  G^(z)  and  G^lz)  be  chosen  as  the  mxm  Wiener-Hopf  factors 
of  the  two  matrix  equations 

G^(z)Gj^(z)  =  l^+H(z)H,;,{z)  ,  (17) 


G2(z)G2.^(z)  =  1^+H,,(2)H(z)  .  (18) 

Clearly,  because  of  the  form  (12)  of  H(z),  both  H(z)H..j(z)  and  H,,j(z)H(z)  are 
paraconjugate  hermitean  matrices  with  powers  of  z  restricted  to  lie  between 
-(n-1)  and  n-1. 
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From  a  classical  result  [3],  G^(z)  and  polynomial  and  uniquely 

determined  up  to  multiplication  on  the  right  by  arbitrary  constant  mxm  unitary 
matrices.  In  addition,  their  degrees  do  not  exceed  n-1  and  they  are  minimiun- 
phase,  i.  e. , 

detG^(z)  •  detG2{z)  ^  0  ,  jz]  <  1  .  (19) 

But,  as  it  is  obvious  that  the  right-hand  sides  of  (17)  and  (18)  are  also 
i  9 

nonsingular  for  all  z  =  e"*  ,  0  real,  we  actually  have  the  stronger  inequality  (13) 
which  Includes  the  boundary  of  the  unit-circle,  jz]  =  1. 


If  we  now  observe  that 


Sf,(p)  = 


G[^z) 

t 

£ 

o 

0 

G‘^z) 

H(z) 

e-^P^l 

m 

e-"P"l 

m 

-z%(z) 

(20) 


it  is  straightforward  to  verify  (16)  with  the  aid  of  (17)  and  (18).  Thus,  Sj^{p)  is 
meromorphic  and  regular  paraconjugate  unitary.^ 

Since  a  right-inverse  is  also  a  left-inverse,  (16)  implies  that 

Sjij^(p)Sj^(p)  =  •  (21) 


In  turn,  (21)  yields' 


-  (G^Gj^)H  =  0^ 


(GjG^,,)-Uh(G2G2*)’^H^=1 


m 


From  (23), 


G'^z)H(z)G2(z)=  L(z)  =  Gj,j,(z)H(z)G'ji(z) 


(22) 

(23) 

(24) 

(25) 


It  follows  from  (13)  that  Gj^^(z)  and  G2^(z)  are  analytic  in  lzl<  1  so  that 
the  domain  of  analyticity  of  Sj^(p)  includes  the  jiu-axis  in  the  p-plane. 

5 

Arguments  are  omitted  wherever  convenient. 
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But  the  right-hand  side  of  (25)  is  analytic  in  |z|_>l  whereas  the  left-hand 
side  is  analytic  in  |z|  <  1  and  this  means  that  I-<(z)  is  analytic  in  the  entire 
finite  z-plane.  By  Liouville's  theorem,  L(z)  is  a  polynomial  matrix  and  the 
proof  of  lemma  1  is  complete,  Q.E.O. 

Comment  1.  If  H(z)  is  chosen  real,  G^(z)  and  G^Cz)  can  also  be  con¬ 
structed  real  and  Sj^(p)  is  a  real  regular  paraunitary  matrix. 

Let  us  now  suppose  that  a  2m-port  N  of  generic  type  defined  in  lemma  1 
is  closed  on  its  output  side  on  a  passive  m-port  load  with  normalized  scattering 
matrix  description  r(z).  Of  course,  r(z)  is  mxm  and  bounded  and  therefore 
analytic  in  |  z  |  <  1 . 

Let  the  four  mxm  partitions  in  (15)  be  labeled  S^^(p),  ^^^(p)  and 


^22^^^*  first  principles  [4], 

S^Sn+Sjjni^-S^^n-'s^l  (26) 

is  the  normalized  scattering  matrix  of  the  resultant  m-port  N__ 

XTj  • 

Thus,  making  the  appropriate  substitutions  and  using  (17), 

S(z)  =  G[^H+z”Gj’^r(l^+z"G‘^H,.r)’^G‘^  (27) 

=  G[^HG2(l^+2*'G'^H.,r)  +  z"r)(G2i-z’^H.j,r)'^  (28) 

=  G'^HG2+  z"(ln,+  HH^)r)(G2+  z%  F)'^  (29) 

=  G'^HG2+z"GjG^,,r)(G2+z"H..r)'^  (30) 

=  (L+z''Gj,.r)(G2+z"H.r)’^  (31) 

=  (L(z)  +a  j(z)  r(z))(G2(z)  +  fi(z)  r(z))  ■  ^  (32) 

where  L(z)  =  G ^  Vz)H(z)G2(z)  and^ 

<3j(z)  =  z‘'g^,^(z)  ,  (33) 


(j^(z)  and  fl(z)  are  the  mxm  matrix  polynomials  reciprocal  to  Gi(z)  and 
H(z),  respectively. 
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(34) 


fl(z)  =  z”h^(z)  . 

Consequently,  since  it  has  been  shown  in  lemma  1  that  L.(z)  is  polynomial, 


s  =  (L+Gjr){G2+flr)"^ 

(35) 

appears  as  a  bilinear  matrix  transform  of  F  with  polynomial  matrix  coef¬ 
ficients  that  depend  on  z  alone. 

Lemma  Z.  Let 

r(z)  =  zp{z) 

(36) 

where  p(z)  is  an  arbitrary  mxm  bounded  matrix  in  z.  ^  Then, 
efficients  in  the  power-series  expansion  of 

the  first  n+1 

S  =  (L+z<3^  P)(G2+  zf?p)‘^ 

(37) 

about  z  =  0  are  independent  of  the  choice  of  p(z). 

Proof,  Let  S  (z)  denote  S(z)  with  p(z)  =0  .It  suffices 
o  m 

to  prove  that 

S(z)  -S^(z)  =  0(z"+^) 

(38) 

in  the  neighborhood  of  z  =  0. 

Evidently, 

S^(z)  =  L(z)G■^z)  =G'^z)H(z) 

(39) 

since  G  j  L  =  HG^.  Thus, 

S-S^  =  (L+z5jP)(G2+  zfip)’^-G'^H 

(40) 

so  that 

G^(S-S^)(G2+  zfip)  =  Gj(L+  zG^p)  -  H(G2+  zfip) 

(41) 

=  z{GjGj-Hfi)p 

(42) 

=  z  p  , 

(43) 

*Such  a  r{z)  is  automatically  bounded.  Moreover,  every  S{z)  generated 
by  (37)  IS  also  bounded  i£  c(z)  is  bounded. 
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from  (17).  Or, 


(44) 


S-S^^z'^+^G'Vc^+zflp)'^ 

However, because  G^(0)  and  G^iO)  are  nonsingular  and  p(z)  is  analytic  in 
lz|<l,  (44)  implies  the  desired  result  (38),  Q.E.D. 

Comment  2.  Since  S(z)  is  bounded, 

Z(z)  =  (l^-hS(z))(l^-S(z))‘^  (45) 


is  positive  whenever  it  exists.  Clearly,  as  is  seen  from  (37)  and  (12), 

S(0)  =  L(0)G‘^0)  =  G‘^0)H(0)  =  0^  (46) 


so  that  neither  1^+ S(z)  nor  l^-S(z)  are  identically  singular. 

Z  ^(z)  are  both  well-defined. 

Thus,  Z(z)  and 

Let  us  write 

n 

Z(z)  =  1^+  2  2C(k)z‘'+  0(z’^'^S  . 
k=l 

(47) 

According  to  lemma  2,  the  mxm  coefficients  C(k),  k=  l-»n,  are  independent  of 
the  particular  choice  of  bounded  p(z)  and  are  therefore  calculable  from  the 

power -series  expansion  of 

(48) 

Since  the  factors  in  (48)  commute  we  can  also  write 

^o  =  <‘m-So''‘lWSo' 

(49) 

(50) 

=  (Gj-H)'^G^+H)  =  2P;^Q^ 

(51) 

where 

P^(z)=G^(z)  .H(z) 

(52) 

and 
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(53) 


Qn(-) 


Gj(2)  +  H(z) 

- z - 


By  direct  calculation  and  Eq.  (17), 

P  (2)Q  (z)  +  Q  (z)P  (z)  =  1 

n  n*'^  n  n^  m 


(54) 


and  we  are  now  ready  for  the  remarkable  lemma  3. 

Lemma  3.  1)  The  polynomial  matrices  Pj^(2i)  and  Qjj(2i)  are  nonsingular 

in  the  closed  unit  circle,  i.  e. , 


det  P^(z)  •  detQ^(z)  ^  0  ,  lzl<l 


(55) 


Z  (2)  +  Z  (z) 

2)  ^  ^  =(P„„(z)PJz)) 


(56) 


3)  The  coefficients  C(k),  k=  l-»n,  determine  the  interpolatory 
2m -port  N  completely  up  to  within  a  frequency  insensitive  all-pass  cascaded  at 
the  output  side. 

Proof.  1)  Suppose  first  that  for  some  z  ,  Iz  I  <  1,  det  P  (z  )  =  0.  Then. 

o  '  o  no 

there  exists  a  nontrivial  m -vector  a  such  that  a^P  (z  )  =  0^  . 

-  -  n'  o'  -m 

Thus,  from  (51)  and  the  analyticity  of  2^{z)  in  Izl  <  1, 

a’i'P^(z^)Z(z^)  =  0;;^=2a«Q^(z^)  .  (57) 

Therefore,  invoking  (52)  and  (53), 

a»(P„(z  J  .  a»G  j(z^)  =  o;_,  (58) 

which  is  impossible  because  G^(z^)  is  nonsingular  in  [zj  <  1. 

On  the  other  hand,  if  [z^j  =  1,  o)  = •  °n#^^o>  =  °J<V 

(54)  evaluated  at  z  =  z  reads 

o 


m-11 


Thus,  by  multiplying  (59)  on  the  left  and  right  by  a*  and  a,  respectively, 
we  obtain  a*a=0,  another  contradiction. 

2)  According  to  (51)  and  (54), 

Z  +  Z^  .  ^  . 

=  Q  P  ^)P~^  =  (P  .  (60) 

2  n  '  n  n*  n  n*!*'  n#  n*  n'  '  ' 


3)  Let 


P  (z)  =  A„+ A,  z  +  .  . .  +  A  z  , 
n'  '  0  1  n  ’ 


where  the  A's  remain  to  be  determined  in  terms  of  the  C's. 
From  (60), 


Zo(^)  +  2^(z) 
- 2 - ^ 


P  (z)  =  z  .P  ^(z) 
n'  n 


where 


P  (z)  =  z”p  j,{z)  =  A*  +  A*  z+...+A*z”  . 
n  '  n*  n  n-1  o 


Part  1)  permits  us  to  conclude  that  Z  (z)  and  P“'^(z)  are  both  analytic 

in  Izl  <  1  and  it  is  therefore  legitimate  to  set  z  =  e'^”,  6  real,  and  to  equate 
'  ”  ikS 

coefficients  of  e*^  on  both  sides  of  (62),  all  k. 

Carrying  this  out  for  k=0-»n  we  obtain  the  linear  system, 

A*XT  =  [1  0  ...  0^]  (64) 

u  n  mm  m 


where 


X  =  [Ao  Ai  ...  A^l 


C(-l)  1 

T  =  m 

n 


C(l)  ...C(n) 

1  ...C(n-l) 


C(-n)  C(-n+l) ...  1 


1 

E 

m 

n 

E* 

T  , 

n 

n-1 

8  * 

For  any  k,  C(-k)  =C’^(k).  Matrices  which  possess  the  structure  in  (66) 
are  said  to  be  normalized  block  Toeplitz. 


Since 


V9) 


Z  (ej®)  +  Z*(ej®) 
o  o  ' 

- z - 


=  (P^(ejV„(ej®))'^ 


(67) 


is  obviously  positive -definite  and  over  -Tr<0<  TT,  and  since  it  also  satisfies 
the  Paley -Wiener  criterion  (5),  it  follows  [6]  that 

T^  =  Tj>0  .  (68) 

The  solution  of  (64)  is  now  routine. 

Let  Aq  be  determined  from  the  Gauss  factorization 

A*A-=(1  -E  T'\e*)‘^  (69) 

0  0m  n  n-1  n'  '  ‘ 


and  let  (T  ^),  denote  the  matrix  formed  with  the  first  m  rows  of  T 

nil 


-1 


Then, 


X  =(A«)‘^  (T;S^.  .  (70) 

Equation  (69)  fixes  Aq  uniquely  up  to  multiplication  on  the  left  by  a 

constant  mxm  unitary  matrix  U.  In  view  of  (70),  Pjj(z)  is  then  also  unique  up 

to  within  the  same  pre-factor.  The  next  step  is  to  produce  a  formula  for  the 

coefficients  of  Q  (z). 

n 

Let 

and  let 

Y  =  [Fq  ...  F^]  .  (71) 

From  (51), 

2Q^(z)  =  P„(z)Z^(z)  (72) 

and  a  comparison  of  coefficients  of  like  power  of  z  ,  k=  0  ->n,  yields 

2Y  =  XZ^  (73) 

where 
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■»**•  -Tt-  • 


2C(1)  ...2C(n) 

1 

m 


Z  = 
n 


.  2C(1) 
1 

m 


(74) 


Thus,  Qj^(z)  is  also  uniquely  determined  up  to  within  the  unitary  premultiplier  U. 

Consequently,  H(z)  =  j  P^(z) -Q^(z)  and  G^(z)  =  ^P^{z) +Q^(z)  are  both 
determined  up  to  multiplication  on  the  left  by  U,  and  then  can  be  found 

as  the  minimum -phase  polynomial  solution  of  (18).  As  mentioned  previously, 
is  unique  up  to  multiplication  on  the  right  by  an  arbitrary  constant  mxm 
unitary  matrix  V. 

If  we  now  return  to  (15),  it  is  immediately  seen  that  Sjj(p)  is  unique, 

Si2(p)  is  unique  up  to  the  right  multiplier  U*,  S2j(p)  up  to  the  left  multiplier 

V*  and  S22(p)  up  to  the  left-right  multipliers  V*  and  U*,  respectively.  But 

this  can  be  interpreted  to  mean  [4]  that  every  compatible  interpolatory  2m -port 

N  is  obtained  from  a  particular  one,N  ,  say,  by  cascading  N  on  the  right  in  a 

P  P 

2m-port  all-pass  with  normalized  scattering  description 


1 

m 

u* 

V* 

1 

m 

This  completes  the  proof,  Q.E.D. 

The  preceding  lemma  suggests  a  very  natural  definition. 


(75) 


Def.  1.  A  set  of  mxm:  matrices  C(k),  k=  l-*n,  is  said  to  admit  the 
strong  interpolation  property  (the  SIP),  if  there  exist  two  mxm  polynomial 
matrices  Pjj(z)  and  Qjj(z)  of  degrees  less  or  equal  to  n  such  that, 

I^.  detP^(z)  .  detQ^(z)  0,  lz|<l; 

^2-  Pn(^)^n*<^)  +  Qn<^)Pn«(")='m  ^ 
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I^.  For  [zj  <  1  , 

n 

Z^(z)  =2P"^z)QJz)=  1^+  Yj  2C(k)z*"+0(z"'^^)  ; 

k=l 

I^.  Z^(z)  is  positive. 

Lemma  4.  A  set  of  mxm  matrices  C(k),  k=  l->n,  admits  the  SIP  iff 


1 

m 

C(l) 

.  .  .  C(n) 

C(-l) 

1 

m 

.  .  .  C(n-l) 

>  0  . 

(76) 

C(-n) 

C(-n+l) 

•  •  •  ^m 

Proof.  Let  the  data  admit  the  SIP.  Then,  from  1^,  I^,  and  the  definition 


of  Z^(z), 


K^(0). 


Z  (ej®)  +  Z*(e->®) 


=  (P*(e-’®)P^(ej®))-^ 


n  n 


(77) 


is  Lj  and  positive -definite  over  -Tr<9<TT. 

Furthermore,  in  view  of  I^* 

^  /^e'j^K^(9)de  =  C(k)  ,  -n<k<n  ,  (78) 

-TT 

and  since  it  is  obvious  that  K^(9)  is  nonsingular  for  all  9,  it  follows  from  the 
discussion  at  the  end  of  Section  II  that  T^>  0. 

Sufficiency,  A  review  of  the  proof  of  lemma  3  shows  that  the  poly¬ 
nomials  Pj^(z)  and  Q^(z)  must  be  given  by  the  formulas, 

P^(z)  =  Xe(2).  QJz)  =  Ye(z)  (79) 


where  X  and  Y  are  determined  by  means  of  Eqs.  (69) -(74)  and 

4  (z)  =  tl__.  zl  .  •  .  z  1  1  •  (80) 

“1  rn  m 
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Of  course,  the  positisre -definiteness  of  guarantees  that  and 

Q  (z)  are  well-defined  and  it  remains  to  be  shown  that  they  possess  all  the 
n 

requisite  properties. 

Note  first  that  I2  is  equivalent  to 


1  X(|  +  Z  i  iJX*  =  1 

2  n  n  m 

But,  as  is  easily  seen, 

e(z)e#(z)  =  l+n(z)  +  f^(z) 


where 


n(z)  = 


m 


zl 


m 


m 


z  1 


m 


m 


zl 


m 


n-1. 


m 


zl  0 
m  m 


(81) 


(82) 


(83) 


Hence,  since 


T  =  4-  (Z  +  Z*  ) 
n  2  n  n ' 


and  JL,  commutes  with  Z  ,  we  find  that 
V  n 

.1. 


■  I  1 


=  T^+^(n+a,)z;+iZn(n+fi,) 


Z  +Z*  Z  +Z"' 
n  n  ,  n  n 


Q 


=  T  +a,T  +  T  n  . 
n  n  n 


(84) 


(85) 

(86) 

(87) 


^To  a  great  extent,  it  is  this  commutativity  that  is  responsible  for  so 
many  of  the  remarkable  properties  exhibited  by  Toeolitz  matrices. 
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*  - 1  -1 

Now,  by  using  the  explicit  formula  X  =  (A*)  •  (T  ),  we  obtain, 

u  n  X , 

in  succession, 

XT^X*  =  (A*)"^  .  (A*)  =  ,  (88) 

XT  n=  (A*)"^.  [1  0  .  .  0  ]n=0  (89) 

n  u  mm  m 

and 

X(T  +a.T  +T  n)X*=l  .  (90) 

Thus,  I2  is  validated. 

The  remciining  three  properties  can  be  disposed  of  almost  simultaneously 
by  proving  that  Z^(z)  =  2P^^(z)Q^(2)  is  positive  in  lz|<  1. 

First,  P^(0)  =  Aq.  But,  in  view  of  (69)  and  the  premise  >  0,  Aq  is 
nonsingular  so  that  detP^(z)^  0.  Thus,  ^^(z)  is  well-defined  and  analytic  in 
some  neighborhood  of  z  =  0. 

By  construction,  its  power  series  expansion  about  z  =  0  must  be  of  the 

form 

00 

Zo(z)  =  1^+  20(10)2^"  (91) 

k=l 

where  C(l),  C(Z),  ...,  C(n)  constitute  the  prescribed  data.  Let  r>  0  denote 
the  radius  of  convergence  of  (91).  Of  course,  since  we  have  not  yet  estab¬ 
lished  that  P^(z)  is  minimum -phase,  i.e.,  that  det  P^(z)  ^  0,  |zl<  1,  it  can¬ 
not  be  asserted  that  r  >  1. 

Nevertheless,  by  equating  coefficients  on  both  sides  of  the  identity 

(  S  A,  )•(*„*  S  2C(k)z'')  =  2  f  F,  s'  (92) 

f=0  k=l  i=0 

we  obtain  the  formulas. 
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^  =  0*^“°  ’  (93) 

where  is  given  by  (74)  with  n  replaced  by  /, 


(The  matrices  and  0^  are  upper-triangular  block  Toeplitz.) 

Our  aim  is  to  prove  that 

+E*)  =  a"\a^  a*)(ap'^  >  o  (96) 


for  all  £  >  0,  Or,  equivalently,  that 

Yjj  s  a|‘  >  0  ,  i  =  0  -»  =0  .  (97) 

By  hypothesis,  >  0  and  (97)  is  correct  for  t  =0-*n.  Clearly  then, 
the  desired  result  >  0,  i>n,  follows  trivially  if  it  can  be  shown  that  for 
all  such  I  , 


^  ^m(f  -n+1)  ^  '^n-1  ‘ 


in-18 


To  accomplish  this  we  shall  exploit  the  already-derived  "even-part" 
condition, 


Pn(z)Qn*(-)  +  =  'm  * 


(99) 


We  leave  it  to  the  reader  to  show  that  a  direct  comparison  of  like 
powers  of  z  on  both  sides  of  (99)  leads  to  the  1-step  update. 


1  • 

n  m  n-i 


(100) 


Since  (100)  is  (98)  with  f  =n,  the  natural  way  to  proceed  is  by  induction. 
Evidently, 


*f +1 


"12 


and  8 


i+1 


F  , 

0 

12 

(101) 


where 


10 


■^12”  +1  ^ 


(102) 


and 


(103) 


By  direct  calculation. 


in  which. 


'^i+l~'^/+l^f +1  ^i+l'^+l  ' 


f +1 

Yu  =  S  <Ai,f5+F^A*) 
k=0 


"ll 

^12 

^2 

(104) 


(105) 


and 


■nr 


For  f  >n.  A-  =  F.  s  0  . 

X  £  ni 


in -19 


If  we  now  take  stock  of  the  previous  footnote,  it  is  quickly  observed 
that  (105)  and  (106)  reduce  to 

n 

+  <107) 

k=0 

and 

.  (108) 

respectively. 

However,  by  equating  the  (1,1)  and  (1,2)  blocks  on  the  two  sides  of 
(100)  it  is  found  that  =  1^  and  ®  whence. 


’'1+1= ‘m* 


(109) 


and  the  induction  step  is  completed.  Thus, 

T^>0,  /=0-»oo  .  (112) 

According  to  a  generalization  of  a  theorem  of  Schur  (Appendix  A), 

(112)  is  a  sufficient  condition  for  the  power-series 

<X5 

ZC{k)z^  (113) 

k=l 

to  define  a  positive  function  of  z  in  ]z|<  1. 

However,  since  this  series  and  Z^(z)  coincide  in  [z]  <  r,  they  must, 
by  analytic  continuation,  coincide  in  all  of  |z|  <  1.  We  have  therefore  shown 
that  Z^(p)  =  2P'^z)Q^(  z)  is  positive  and  it  remains  to  check  property  Ij. 
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w 


The  positivity  of  Z^{z)  automatically  guarantees  its  2uialyticity  in 

jzj  <  1.  Hence,  as  shown  in  lemma  3,  if  for  some  z^,  jz^l  <  1,  there 

exists  a  nontrivial  vector  a  such  that  a*P  (z^)  =  0^  ,  then  necessarily, 

a*Q  (z  )  =  o'  . 

»  n  o  •m 

Consequently,  from  (99), 


a*P  (z^)Q  (z  )  +  a«Q  (zJP  (z  )  =  =  0 ' 

••  n  o  n  o  ->  n  o  n  o  o—  —m 


(114) 


which  implies  that  z  =0.  But  this  is  a  contradiction  because  P  (0)  =  A  is 

o  n'  '  o 

nonsingular. 

Similarly,  since  it  follows  immediately  from  Eq.  (59)  that  Pjj(®)  and 
Q^(z)  are  both  nonsingular  on  the  bovmdary  of  the  unit  circle,  the  assumption 
jz^l  =  1  is  equally  untenable.  Finally,  by  observing  that 

z;^(z)  =  iQ;\z)p^(z)  (115) 

is  also  positive,  we  can  apply  exactly  the  same  reasoning  to  establish  the 
non  singularity  of  Q^(z)  in  1,  Q.E.  D. 

Corollary.  Let  the  mxm  matrices  C(k),  k=  l-*n,  admit  the  SIP,  Then, 
there  exists  an  interpolatory  2m-portN  such  that 

n 

Z^(z)  =  1^+  2C(k)z‘'+0(z“+^)  .  (116) 

k=l 

Furthermore,  N  is  unique  up  to  a  frequency-insensitive  2m-port  all-pass 
cascaded  on  its  output-side. 

Proof.  Let  us  employ  the  procedure  des  r-i’- -  d  in  lemma  4  to  construct 

P  (z)  and  Q  (z)  from  the  prescribed  data.  The  'i  know, 

n  n 

det  P^(z)  .  detQ^(z)  ^  0  ,  lzl<l  ,  (HI) 


and  Z  (z)  s  2P  ^ 
o'  n 


(z)Q^(z)  is  positive. 
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As  before,  let 


G^(z)  =  ip^(z) +  Q^(z)  (118) 

and 

H(z)  =  i  Pn(2)  -  Qn(2)  .  (119) 

Since  property  is  satisfied, 

G^(z)G^,^(z)  =  l^  +  H(z)H^(z)  (120) 

and  it  remains  to  prove  that  G^(z)  is  minimum -phase. 

Obviously,  the  positivity  of  Z^(z)  implies  that  of  +  and  this 

entails  that  the  latter  be  nonsingular  in  |z|  <  1.  But  in  view  of  (118), 

detGj(z) 

det(l^  +  2^(z))=2-.g^^^  (121) 

and  it  follows  that  detGj(z)  is  free  of  zeros  in  jzj  <  1. 

The  construction  of  N  is  now  completed  by  determining  02(2)  as  the 
minimum -phase  solution  of  (18).  Of  course,  the  arbitrary  constant  unitary 
matrix  right -multiplier  V  can  again  be  absorbed  in  an  output  2m-port  all¬ 
pass,  Q.  £.  D. 

Comment  3.  The  recursion  in  (98)  was  discovered  (in  the  scalar 
case)  by  the  author  and  published  for  the  first  time  in  Reference  7.  It  is 
another  manifestation  of  the  incredibly  rich  substructure  exhibited  byblock- 
Toeplitz  matrices. 

It  is  also  important  to  understand  that  the  integer  n  is  det>  -mined 

by  the  number  of  data  and  serves  only  as  an  upper  bound  on  the  degrees  of 

P  (z)  and  Q  (z).  For  example,  the  set  C(k)  =  0  ,  k  =  l-»n,  admits  the  SIP 

n  n  m 

since 

T^=1>0  .  (122) 

'■^Actually,  detC|(z)/0in  1*|<  1. 
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Clearly,  in  this  case 


and 


(T‘h,  =  Cl  0  .  .  .  0  1 

'  n  '1 ;  mm  m 


(123) 


£ 

n 


[0^  0^  ...  0  ] 
mm  m 


(124) 


so  that 


A 


0 


=  1 

m 


(125) 


and 

P^(z)  =  I  =  2Q  (2)  (126) 

n  m  n 

are  compatible  choices.  Thus,  Pjj(z)  and  0^(2)  are  constants  and  total  de¬ 
gree  reduction  has  occurred. 

At  this  stage  we  are  ready  to  state  and  prove  our  first  major  theorem 
which  serves  as  the  foundation  stone  for  the  multichannel  spectral  estimation 
technique  developed  in  the  next  section.  However,  we  must  first  insert 
another  definition. 

Def.  2.  A  matrix  Z(2)  is  said  to  be  full-passive  if  it  is  positive 

and 

K(9,  .  I  (127) 

is  nonsingular  for  almost  all  9  in  -TT <  6<  TT ;  i. e. ,  if 

detK(9)^0,  a.e.  (128) 

Similarly,  an  mxm  matrix  T(z)  is  full-bounded  if  it  is  bounded  and 

det(l^-  r*(ej®)r(e^®)  0  ,  a.e.  (129) 

Theorem  1.  Let  C(k),  k  =  0  -♦n,  be  any  given  set  of  n+1  constant 

mxm  matrices  such  that  C(0)  =  1  and  T  >  0,  and  let  N  denote  tlie  (unique) 

m  n 

associated  interpolatory  2m -port  defined  in  the  corollary  to  lemma  4. 


Ill- 2.3 


Let  P  stand  for  the  collection  of  all  mxm  full-passive  matrices 
n 

n 

Z(z)  =  1^  +  S  2C(k)z‘^+  .  (130) 

k=l 

(Any  such  Z(z)  is  said  to  interpolate  to  the  prescribed  data.) 

Then,  the  class  ^  coincides  with  the  set  of  all  mxm  positive  matrices 
Z(z)  generated  by  formulas  (45)  and  (37)  in  which  p(z)  is  an  arbitrary  mxm 
full -bounded  matrix. 

Proof.  For  n=0,  is  the  set  of  all  full-passive  mxm  matrices 
Z(z)  such  that  Z(0)  =  1^.  Hence,  in  particular, 

S(z)=(Z(z)  -  l^)(Z(z)  +  1^)"^  (131) 

12 

is  bounded  and  equal  to  0^  for  z  =  0.  Thus,  by  Schwartz’s  lemma, 

S(z)  =  zp{z)  ,  (132) 

p(z)  boxinded. 

In  other  words,  all  positive  Z's  that  satisfy  the  interpolation  con¬ 
straint  Z(0)=  1^  are  obtained  by  letting  p(z)  range  over  all  bounded  mxm 
matrices  in  the  formula, 

Z=  (l^  +  zp)(l^-zp)'^  .  (133) 


But  Eq.  (132)  is  (37)  with  H(z)  =  0^  and  G^(z)  =0^(2)  =  1^  so  that  the 
corresponding  interpolatory  2m -port  N  has  the  normalized  scattering  description 


- 


■  0 

1 

m 

m 

1 

m 

0 

m 

(134) 


Algebraic  manipulation  of  Eq.  (133)  yields, 

(1  -ej®p(ej®))*K(e)(l  -ej®p(e^®))  =  I  -  p*(ej®)p(ej®)  .  (135) 


Let  r(2)  be  an  mxm  bounded  matrix  such  that  r(0)=  0^. 
r(z)/z  is  also  bounded. 


Then  [9], 


Consequently,  Z(z)  is  full-passive  iff  p(z)  is  full-bounded  and  the  case 
n=0  is  verified.  For  n>l  we  proceed  by  induction. 

12' 

First,  Z(z)  obviously  implies  that  Z(z)  Let  the  2m-port 

interpolatory  cascade  N.  constructed  with  the  data  C(k),  k=  0  -»n-l,  be 

^  1 
parametrized  as  described  in  lemma  1  by  the  mxm  polynomial  matrix  H^(z). 

Of  course,  degree  H^(z)<n-1  and  ^11^*^  ^12^*^  denote 

the  companion  left-right  mxm  polynomial  Wiener-Hopf  factors: 


^11^11*  "  ^m'*'  ^1*^1*  ' 


(136) 


G,,G  =1  + 

12  IZ’!'  m  1»  1 


(137) 


Also,  let 


Lj(z)  =  g;;(z)Hj(z)Gj^(z) 


(138) 


The  induction  hypothesis  permits  us  to  assert  that 


where  (Eq.  (3?)), 


Z  =  (1  +  S)(l  -  S) 

m  m 


S  =  (L^+zGjjTi)(Gj2  +  2l5i'n)’^ 


(139) 


(140) 


and  T|(2)  is  some  full-bounded  mxm  matrix.  Thus,  making  the  appropriate 
substitutions  in  Eq.  (44)  and  replacing  n  by  n-1  we  obtain. 


in  which 


S(.)-S  (z) 


Sio(z)  =  G';(z)Hj(z) 


(141) 


(142) 


According  to  lemma  2,  the  power-series  expansions  of  S(z)  and  S^^(z) 
about  z  =  0  agree  up  to  z*^'^  and  this  means  that  the  quotient  on  the  left-hand 

- TT - 

Tn  >  0  implies  that  T^.i  >  0.  Incidentally,  we  already  know  from 
lemma  4  that  empty. 


in  (8). 


This  Hj(z)  must  not  be  confused  with  the  constant  matrix  Hj  appearing 
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I 


side  of  (141)  is  analytic  in  |z|  <  1.  However,  (139)  determines  the  power - 
series  expansion  of  S(z)  up  to  z*^  because  that  of  Z(z)  is  known  up  to  z^. 

If  we  couple  this  with  the  observation  that  S^^(z)  is  completely  fixed 
by  the  data  C(k),  k=0  -»n-l,  it  is  concluded  from  (141)  that 


T1(0)  =  G^^(O) 


S(z)-S^^(z) 


n 


Gi,(0) 


z=0 


(143) 


is  a  known  quantity.  Hence,  ri(z)  must  be  selected  from  the  class  of  mxm 
bounded  matrices  which  at  z  =  0  interpolate  to  the  prescribed  value  given 
in  (143). 

Clearly,  ■n(z)  bounded  implies  that  ll'n(0)ll  <  1  but  actually,  1111(0)11  <  1. 
Suppose  otherwise.  Then,  1^  -  11*(0)T1(0)  is  singular  and  there  exists  a 
vector  a  of  unit  norm  such  that 

a  =  t1*(0)n(0)a  .  (144) 


Or, 


a  =  Tl*(0)b 


(145) 


where 


b  =  Ti(0)a 


(146) 


Evidently, 

a*a  =  b*b  =  1  . 


(147) 


Construct  two  constant  mxm  unitary  matrices  V  and  U  which  incorporate 
a  and  b,  respectively,  into  their  first  columns  and  let 


6(z)  =  U*Ti(z)V  . 


(148) 


The  matrix  5(z)  is  mxm  and  bounded  and  6^  ^(0)  =  1. 

Since  all  elements  of  6(z)  are  analytic  and  bounded  by  unity  in  jzj  <  1, 
we  conclude  by  maximum -modulus  that  6^  ^(z)  =1.  In  turn,  the  requirement 

l^-6#(z)6(z)  >  0  ,  lzl<  1  ,  (149) 
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then  forces  all  the  other  entries  in  the  first  row  and  column  of  6(z)  to  be 
identically  zero  so  that 

6(z)  =  1^  +  6j(z)  ,  (150) 

6j(z)  an  (m-l)x(m-l)  bounded  matrix.  Thus, 

=V*(1^-6*6)V  =  V*(0^+  (1^.1  -6*5^))V  (151) 

is  clearly  singular  for  all  lz|  <1  which  contradicts  the  full -bounded  character 
of  ri(z).  Consequently, 

lj„-‘n=*‘(0)Tl(0)  >  0  (152) 

and  the  renormalization  scheme  described  in  Appendix  B  becomes  available. 
Let  and  R^  be  any  pair  of  hermitean  matrix  solutions  of  the  Eqs. , 

=  i„-n(0)n*(0)  ,  (153) 

R^  =  -n*(o)n(o)  .  (154) 

Then,  from  Schwartz's  lemma  and  the  renormalization  theorem  in  Appendix  B, 

Rj^{lm-’^^*(0»’^(‘n-Tl(0))R^^  =  zp(z)  ,  (155) 

p(z)  an  mxm  bounded  matrix.  In  addition,  it  is  seen  from  Eq.  (BIO)  that  'n(z) 
is  full-bounded  iff  p(z)  if  full -bounded.  The  last  step  is  to  solve  (155)  for  11(2), 
substitute  into  (140)  and  then  rearrange  properly. 

The  end  result  is 

S  =  (L  +  zGjP)(G2+zHp)’^  (156) 

where 

H=  r2^H^+zti(0)Gj2)  =  R;^Gj^  +  zn(0)L:^)  ,  (157) 
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’I 


L=  (Lj+zGjjTi(0))R^^  ,  G2  =  (Gj2  +  zHiTi(0))R;^  .  (158) 

The  proof  that  G^(z)  and  C^(x)  are  minimuin -phase  is  easy.  For  ex¬ 
ample,  since 

Gi  =  R’^l^  +  zn(0)LiG‘^j)Gii  ,  (159) 

it  is  clear  that  G,  is  minimum-phase  iff  the  middle  factor  in  (159)  is  minimxim- 

^  -1  *"  -1 
phase.  But  Lj  =  G  j  j^H^G  j2#  so  that  =  G ^2  H ^G ^  ^  and 


Ll°u=<=u«l 


(160) 

»  1  ^ 

I  /  _  VT  T 


However,  from  (137)  and  maximum -modulus  it  is  seen  that  G^2{2)^l(2) 
is  bounded.  Thus,  since  ||zr|(0)||  <  1,  for  |z|  <1, 


-1. 


1^+  zn(0)L^(z)Gj;(z) 


(161) 


is  a  positive  matrix  and  its  determinant  is  necessarily  devoid  of  zeros  in 
|z|<  1.  A  similar  proof  applies  of  course  to 

Equations  (157)  and  (158)  can  be  combined  into  the  single  matrix  version, 


^1 

H 

^z* 

=  T(z)  . 


^11 

«1 

^1* 

*^12# 

(162) 


where 


T(z)  =  (R^^  +  R^^  . 


1 

m 

z"t1(0) 

2’”l1*(0) 

1 

m 

(163) 


Let 


1 

0 

m 

TD 

0 

-1 

m 

m 

Then,  by  direct  verification, 
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(164) 


and 


T  OpTj^  -  CJj 


^1*^1  '  ^12*^12  ~  2*^1**^  11*^  11*^1*^  12*  ”^12**^12 

=  Gi2:i.(Hj*Hj-Gj2Gj2^)G-J*=  -1^  .  (164a) 


Hence, 


^1  ! 

H 

/T 

■°1 

H 

^  1 

^2* 

^2* 

=  o. 


J* 


(165) 


1.  e. , 


and 


^2^2*  = 


L  =  G-^HG2  =  G^^HG-i  . 


(166) 

(167) 

(168) 


In  short,  it  has  been  shown  that  every  member  of  ^  can  be  synthesized 

n 

by  an  appropriate  full-bounded  termination  of  the  output-side  of  the  2m -port 
1 4 

N  defined  by  H(z).  Since  part  3)  of  lemma  3  assures  us  that  this  2m -port  is 
an  essentially  unique  construct  from  the  full  set  of  coefficients  C(0),  C(l),  ... 
the  proof  of  theorem  1  is  complete,  Q.E.D. 

In  the  next  section  the  procedure  developed  in  theorem  1  will  be  given  a 
condensed  algorithmic  formulation  and  then  applied  to  the  problem  of  spectral 
estimation. 


T? 


Note  that  degree  H(z)  £n  and  H(0)  = 


0  . 

m 


C(n). 
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IV.  THE  CLASS  OF  ALL  SPECTRAL  ESTIMATORS 

Let  us  now  suppose  that  we  have  available  n+1  error -free  mxm  co- 
variance  samples  C(0),  C(l),  C(n)  of  an  m -dimensional  full-rank 
random  vector  process  x^.  We  shall  produce  a  parametric  formula  for  the 
class  of  all  spectral  densities  K(6)  consistent  with  the  given  data. 

Theorem  2  (The  Master  Result).  Let  x^  denote  a  full-rank  m-dimensional 
random  vector -process  and  suppose  that  its  first  n+1  matrix  mxm  covariance 
samples  C(0),  C(l),  . . . ,  C(n)  are  known  without  error. 


Let  T^  and  T^(-l)  denote  the  Toeplitz  matrices  constructed  on  C(0), 

C(l),  C(n)  and  on  the  "reversed"  data  C(0),  C{-1),  C(-n),  respectively. 

(Clearly,  T^(-l)  is  obtained  by  replacing  each  block  C(k)  in  T^  by  C(-k),  k=  0-»n.) 
Let 


|(z)  =  (1  zl 
*  m  m 


n  .  / 
=  1  ) 


(169) 


Then,  the  following  3-step  algorithm  generates  all  spectral  density  matrices 
K(G)  compatible  with  the  given  data. 

1)  Affect  the  (unique)  Gauss  factorizations, 

T  =A*A  ,  T  (-1)  =  A*(-1)A  (-1)  (170) 

n  n  n  n  '  n  '  n  '  ' 


where  A^  and  A^(-l)  are  both  lower -triangular  with  positive  diagonal  scalar 

entries.  Let  M  and  M  (-1)  represent  the  two  submatrices  formed  with  the 

- 1  -1 

first  m  columns  of  A^  and  A^  (-1).  respectively. 

2)  Define  two  mxm  polynomial  matrices 

P^(z;l)  =  M*  |(z)  (171) 

and 

PJz:-!)  =  |'(z)M^(-l)  ,  (172) 


of  degrees  <  n  and  let 
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^^P^(z;-1)  =  z"p^,,(z:-l). 
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(173) 


D^(z)  s  P^(z;l)  -  zp(z)P^(z;-l)  . 


3)  Then, 

K(0)  =  D;^e^®)(l^-p(ej®)ff!‘(eJ®))D;^(e-'®) 


(174) 


where  p(z)  is  an  arbitrary  full-bounded  mxm  function. 

Proof.  Assume  that  C(0)  =  1  and  let  us  show  that  P  (z;l)  =  P  (z). 

m  n  n 

As  we  known,  Pjj(z)  =  X£(z)  where 

X  =  (A*)‘^  .  (T'S^.  (175) 

and 


AtA„  =  (1  -E  T‘\e*) 

00  ni  nn-ln 


,-l 


(176) 


But  clearly,  from  T*^  =  A^^(A’>^)*^  and  the  lower-triangular  character 


of  A^,  it  follows  immediately  that 


(T‘^),  =A:i.M* 

n  '1;  00  n 


(177) 


where  Aqq  is  the  upper  left-haud  corner  mxm  piece  in  A^.  Furthermore, 
since  the  matrix  on  the  right-hand  side  of  (176)  is  the  upper  left-hand  corner 
mxm  piece  in  ,  we  can  also  make  the  identification. 


Thus, 


(178) 


M*  =  M* 
n  n 


(179) 


and  r’jj(z)  =  M*£(^)  =  Pj^(z;l).  To  arrive  at  an  understanding  of  the  meaning  of 
P^(z;-1)  we  must  derive  a  general  expression  for  K(0). 

Let  C(k)  denote  the  covariance  function  of  the  full-rank  process  x^.  Then 
necessarily, 

- ^ -  ifl 

D  .(eJ®)  =  I^(eJ®),  0  real. 


(180) 


T.  >  0  and  T.  (-1)  >  0  , 

K  K 

k  =  0-»<»  . 

Hence,  invoking  the  multivariable  Schur  theorem, 

oo 

Z(z)  =  1^  +  S  2C(k)z^ 
k=l 

is  a  positive  mxm  matrix  . 

Clearly, 

K(8).f  C(k)ej“  = 

k=-oo 

where 

Z{e^^)  s  limit  2(re^®)  . 
r'+l-O 

To  describe  all  spectral  densities  K(0)  consistent  with  the  data  it  is 
sufficient  to  exhibit  a  formula  for  the  hermitean  parts  of  all  full-passive 
mxm  matrices, 

00 

Z(z)  =  1^  +  S  2C(k)z*"  +  0(z"'^S  .  (184) 

k=l 

As  a  preliminary  step  we  shall  use  lemma  Z  to  generate  all  corresponding  mxm 
bounded  functions, 

S(z)  =  (Z(z)  -  l^)(Z(z)  +  1^)'^  .  (185) 

Indeed,  every  such  S(z)  is  given  by  Eq.  (37)  with  p(z)  an  arbitrary  full- 
bounded  mxm  matrix.  Nevertheless,  it  turns  out  to  be  more  convenient  to 
work  with  the  equally  valid  expression 

S(z)  =  (Gj+  zpL)‘^(H  +  zpG^)  *  (186) 

that  is  obtained  by  reducing 


(181) 


(182) 


(183) 
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(187) 


instead  of  (26). 

A  straightforward  calculation  with  the  aid  of  (186)  yields,  for  [z]  =  1, 
(G^+  zpL)(l^ -SS*)(G^  +  zpL)*  =  1^ .  (188) 

Or,  since 

^4^  =  (l^-S)‘\l^-SS*)(l^-S*)'^  (189) 

and 

l^-S  =  (Gi  +  zpL)‘^(Gi-H-zp{G2-L))  ,  (190) 

we  obtain 

=  ,  |*i  =  l  ,  (191) 

where 

D^(z)  =  G^(z)  -  H(z)  -  zp(z)(G2(z)  -  L(z))  .  (192) 

Hence, 

D^(z)  =  P^(z;l)  -  zp(z)PJz;-l)  (193) 

if  it  can  be  demonstrated  that 

P^(z;-1)  hG2(z)  -  L(z)  (194) 

agrees  with  (172). 

From  Eqs.  (25)  and  (39), 

LG'*  =  G'*H  =  (195) 

so  that 
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Thus, 


=  •  O’*'' 

z  +  z 

-S^-°*  =  ((G^  -  L)(G^^  -  L^))-^  ;  (197) 

i.  e. , 

.  ^  Z  +Z  ^ 

(G^-L)  .  9.^  .^=  z^.  (G^-L)~^  .  (198) 

Since 

G^^{z)G^{z)  =  l^  +  L^{z)L{z)  .  (199) 

it  follows  that  llS^(2)||  <  1,  lzj<  1.  Hence, 

(1^  -  S^(z))G^{z)  =  G^{z)  .  L(z)  (200) 

is  nonsinguiar  in  (z  (  <  1  and  the  exact  same  reasoning  used  in  lemma  3  to  derive 
(171)  also  yields  (172)  when  applied  to  (198), 

To  complete  the  proof  we  must  remove  the  restriction  C(0)  =  1^  and 
this  is  accomplished  by  exploiting  the  following  series  of  observations. 

O^.  There  exists  a  lower-triangular  matrix  with  positive  diagonal 
entries  such  that 

(R*)'^C(0)R;^  =  1^  .  (201) 

O^.  The  use  of  (174)  in  conjunction  with  the  normalized  data 

(R*)‘^C(k)R"\  k=0-»n  ,  (202) 

is  tantamount  to  reconstructing  all  admissible  normalized  spectral  densities, 

(Rj)'^K(e)R"^  .  *  (203) 


III-34 


O^.  The  lower-triangular  factorizations  associated  with  the  "forward" 
and  "reverse-direction"  Toeplitz  matrices  constructed  on  these  normalized 
samples  equal,  respectively,  ^^^(1)  and  A^(-l)  multiplied  on  the  right  by 


Thus, 


^o^  =  <^o+ 


M  (1)  -♦  A  M  (1)  , 
n  '  on' 

M  (-1)  -♦AM  (-1)  . 
n  '  on''* 


P^(z;l)-  M#(1)A*  |{z)  . 
P^(z;-1)-  |'(z)A^M^(-l) 


(204) 


(205) 


(206) 


(207) 


(208) 


DJz)-»M*(l)A*£(z)-zp(z)Mj;(-l)Aj|'(z)  h  D^(z)  .  (209) 


O..  Since 
4 


we  see  that 


(Rp-W  =  D;‘(i,^-pp»)D,;i 


(210) 


(211) 


where 


(212) 


Q.E.D. 


Comment  4.  The  simplicity  of  the  general  formula  (174)  is  truly  re¬ 
markable,  especially  when  it  is  realized  that  P^(z;l)  and  P^(z:-1)  are  inti¬ 
mately  related  to  the  right  and  left  matrix  orthogonal  polynomials  defined  by 
the  sequence  C(k),  k=  0  -♦n.  This  means  [6,  10]  that  both  can  be  generated 
recursively  (and  simultaneously)  by  readily  available  and  highly  efficient 
algorithms  of  the  Levinson  type,  thereby  obviating  any  need  to  compute  the 
inverses  a'^  and  a‘^(-1). 


Although  we  intend  to  resume  our  studies  of  these  matters  in  depth 
we  shall  nevertheless  give  the  explicit  connection  formulas  for 
the  benefit  of  the  reader  who  wishes  to  proceed  directly  to  the  numerical 
implementation  of  (174), 

Namely,  if  and  denote  the  upper  left-hand  corner  mxm 

blocks  in  and  A^{-1),  respectively,  and  if  Q^”^(z)  and  R^”^(z)  represent, 
in  the  same  order,  the  monic  left  and  right  degree-n  matrix  orthogonal 
polynomials  alluded  to  above,  then  [6  3, 

Q<"Nz)  =  (212a) 

and 

r("^(z)  =  a|j”>*(-1)P^(z;-1)  .  (212b) 

Consequently,  Eqs.  (136)-(142)  in  Reference  6  become  immediately  available 
after  a  suitable  1-1  map  of  some  of  the  symbols. 
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V.  THE  MULTICHANNEL  FLAT-ECHO  ESTIMATOR 


Under  the  assumption  that  is  a  full-rank  random  vector-process  with 
absolutely  continuous  spectral  function,  a  knowledge  of  its  covariance  samples 
C(0),  C(l),  ....  C(n)  permits  us  to  assert  that  admissible  interpolatory 
K(0)  are  given  by  (174)  where  p(z)  is  any  full-bounded  mxm  function. 

Evidently,  as  is  seen  from  this  expression,  K(0)  satisfies  Paley-Wiener 
iff  the  same  is  true  of 

lm‘ • 

This  grcinted,  there  exists  [6]  an  mxm  bounded  matrix  function  Bp(2)  which  is 
analytic  together  with  its  inverse  in  jz]  <  1,  such  that 

1^-  p(e-'®)p*(e-’®)  =  Bp(e-»®)B*(e-’®)  (213) 

for  almost  all  0  in  -tr^  0<  TT, 

From  (213)  and  (174)  it  follows  that 

K(0)  =  B(ej®)B*(e-’®)  (214) 

where 

B(2)  =  D;^z)Bp(z)  (215) 

is  also  analytic  together  with  its  inverse  in  |z|  <  1.  Hence,  formula  (215) 
constitutes  a  complete  parametrization  of  the  set  of  all  admissible  spectral 
density  Wiener -Hopf  factors  consistent  with  the  given  data,  in  terms  of  the 
free  bounded  matrix-function  parameter  Bp(z). 

Since 

B(0)  =  D'^O)Bp(O)  =  Pj^(0;l)Bp(0)  =  P^(0)Bp(0)  (216) 

and 

1  -B„(0)B*(0)  >  0  , 

tn  p'  '  p'  '  —  * 

— n - 

Remember  that  P„(z)  and  P  (z;l)  are  identical. 


(217) 


we  see  that 


P^<0)Pj(0)  -B(0)B*(0)  =  P^(0)(l^-Bp(0)Bj(0))P^(0)  >  0  .  (218) 


In  particular. 
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ldetP^(0)l>  |detB(0)l  .  (219) 

Furthermore,  because  P  (0)  is  nonsingular,  equality  in  (219)  is  possible 
iff  Bp(0)  is  a  unitary  matrix.  But  then,  it  follows  routinely  from  maximum- 
modulus  that  Bp(z)  =  Bp(0)  for  all  z  and  this  in  turn  implies  that 

p(z)  s  0^  .  (220) 

According  to  (173)  and  (174),  p(z)  £  0^  corresponds  to  the  choice  of 


estimator 


(221) 


J9 


and  invoking  the  result 

1  ^ 

2inldetP^(0)| +2<n|detBp(0)j  =2ln|detB(0)|  =  ^  J Jn detK( e)de  £Entropy(K)  , 


-TT 


(222) 


it  is  immediately  concluded  that  Kj^g(0)  is  the  maximum-entropy  estimator 
[11  ].  In  other  words,  obtained  by  the  termination  of  the  basic 

interpolatory  2m -port  N  in  a  matched  m-port  load,  r(z)  =  0^. 

Alternatively,  it  is  seen  from  (174)  that  the  entropy  of  any  admissible 
K(0)  associated  with  the  choice  of  bounded  function  p(z)  is  given  by,^® 

TT 

Entropy(K)  =  Entropy(Kj^2)  '  ^  p*(ej®))de<  Entropy  (K^^.). 

-TT 

(223) 

— W, 
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that 


if  A  and  B  are  both  hermitean  nonnegative -definite,  det(A+ B)  >det  A+ detB. 
Apply  Cauchy's  theorem  to  the  analytic  function  indetB(z)  and  observe 

detK(0)=  ldetB(ej®)l^  . 


^%ince  1,„- p(e'^®)  p*(e^®) ^  0,  a.e.,  the  reciprocal  of  its  determinant 
number  which  is  greater  than  unity  for  almost  all  6. 


IS  a 
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jg 

Clearly,  equality  is  possible  iff  p(e‘^  )  =  0  a.e,  in  i.e,,  iff 

K(9)  =  K^,J.(6). 

In  view  of  (223),  it  is  legitimate  to  define 

^indet'^l^-p(ej®)p*(ej®))  (224) 

to  be  the  entropy  echo-loss  density  at  angle  6.  A  flat-echo  estimator  (an 
FEE)  is  one  for  which  this  density  is  a  constant  independent  of  0. 

One  particularly  simple  (and  numerically  effective  way)  to  generate 
FEE'S  is  to  choose 

p(z)  =  iJid(z)  (225) 

where  |i  is  an  arbitrary  constant  mxm  matrix  such  that 

1  -m/>  0  .  (226) 

m 

and  d(z)  is  an  arbitrary  rational  mxm  regular  paraconjugate  unitrary  matrix. 
Thus,  for  this  subset  of  FEE'S, 

Kj.j;(e)  =  D'^ej®)(l^-WU'’‘)D;JJ(e-'^)  (227) 

where 

D^(z)  =  P  (z).  zud(z)P  (z:-l)  ,  (228) 

n  n  n 

and 

Entropy(Kp^)  =  Entropy(Kj^j.)  -fndet  ^(l^-UH*)  .  (229) 

Note,  that  the  entropy  of  Kpj,(0)  can  be  fixed  in  advance  by  an  appropriate 
assignment  of  M  and  d(z)  can  then  be  designed  to  accomplish  the  necessary 
tuning.  The  latter  emerges  as  a  problem  of  interpolation  with  regular  mxm 
matrix  all-passes. 

Although  a  complete  (recursive)  solution  has  already  been  obtained  for 
the  case  m  =  1  [l  ],  much  work  remains  to  be  done  for  general  m. 


IT 


d(z)d,j,(z)  =  d(i(z)d(z)  =  Ijn 


and  d(z)  is  analytic  in  |  z  |  ^  1 . 
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APPENDIX  A 


The  Multivariable  Schur  Theorem.  The  power  series 


00 

Z(z)  =  +  S  2C{k)z*'  (Al) 

k=l 

defines  an  mxm  positive  matrix  in  |zl<  1,  iff 

Tj^>  0  ,  k=  l-><«  .  {A2) 

Proof.  Suppose  that  Z(z)  is  positive.  Then,  by  definition,  Z(z)  is 
analytic  in  ]  z  |  <  1  and 


Z(z)+Z*(z) 

2 


|zl  <  1  . 


(A3) 


Hence,  for  every  fixed  r,  0<  r<  1,  and  for  every  mxm  trigonometric 
polynomial  matrix 

f(e)  =  fQ+fje"j®+... +fj^e‘j‘'®  ,  (A4) 

9  real,  k>l, 

J^f*(9)  »  .  f(9)de>0  .  (A5) 

-TT 


Or,  carrying  out  the  integration, 

k  .  . 

J  rl^ .  £5C(q-i)fq>  0  (A6) 

i  ,q=0 

where  C(0)  =  1^. 

In  the  limit  as  r  -» 1-0  with  k  held  fixed,  (A6)  goes  into 
k 

J  f^C(q.i)fq>0  (A7) 

i  ,q=0 


and  because  of  the  arbitrary  character  of  the  mxm  matrices  fQ>^^i  *  *  *  • 
we  conclude  that  0. 

Sufficiency.  From  T^>0  we  deduce,  in  particular,  that 


1 

m 

C(k) 

C«(k) 

1 

m 

i.e.,  for  every  k>0, 

llC(k)ll  <1  . 


(A8) 


(A9) 


Thus,  if  Z(z)  is  defined  by  (Al)  we  obtain  the  estimate. 


llZ(z)ll<  1  +  ^  2llC(k)ll  .  Izl*"  <  (AlO) 

k=l 


<  1+  f  Zlzl*" 

k=l 


1  +  I  z  1 
1  -  jz 


<  00  , 


(All) 


|z|  <  1.  Consequently,  Z(z)  is  analytic  in  |z|  <  1  and  it  only  remains  to 
demonstrate  its  positivity. 

This  can  be  achieved  by  observing  [5]  that 

=  y  ( y  c{k)z*^+  y  c(.k)z^)  (aiz) 


o 

it 

k=0 

k=l 

00 

00 

00  00 

s 

i=0 

k=0 

f=0k=l  (A13) 

00 

00 

00  00 

s 

J  C(a-e)z®z® 

+  1  V  C(a-S)z*l9 

0=0 

a=p 

aao  Ma+l 

(A14) 
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.<n- 


where 


=  2  C(a-0)a®z^ 

a,  p=o 


k 


=  limit 
k-»00 


2  l|C(a.e)f^ 
a,  p=0 


>0 

(A15) 


£ 


a 


a 

z 


1 

m 


» 


a  =  0 -»oo  . 


(A16) 


(In  reaching  (A15)  we  have  used  (A7)  and  have  also  noted  that  the  first 
sum  in  (A14)  contains  all  the  elements  on  and  above  the  main  diagonal  of  the 
doubly -infinite  array 

{C(a-e)z'^r^}  ,  a,  0  =  0- 00  ,  (A17) 


while  the  second  contains  all  those  strictly  below.) 
Since  1  -  |zl^>  0  for  [zj  <  1,  (A15)  yields 


Z(z)-t-Z*(z) 

2  _  ^ 


Izl  <  1 


(A18) 


amd  Z(z)  is  positive,  Q.E.O. 
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APPENDIX  B 


Our  proof  of  theorem  1  made  use  of  the  following  interesting  result 
long  known  to  network  theorists. 


The  Renormalization  Theorem.  Let  be  any  constant  mxm  matrix 
such  that  llC^II  <  1.  Let  and  denote,  respectively,  any  pair  of 
hermitean  mxm  matrix  solutions  of  the  two  equations. 


and 


a 


1  -C  C* 
m  o  o 


1  -C'C 
m  o  o 


(Bl) 


(B2) 


Then,  if  W(2)  is  an  mxm  bounded  matrix,  so  is 

r(z)  =  R  (1  -WC*)‘^W.C  )R'^ 

'  a  m  o'  '  o  b 


Proof.  The  obvious  identity 

C*(l  -C  C*)'^ 
o'  m  o  o 


[1  -C*C  )‘^C* 
'  m  o  o  o 


yields. 


C*R"^  =  RT^C* 

o  a  bo 


(B3) 


(B4) 


(B5) 


Hence,  upon  multiplying  both  sides  of  (B5)  on  the  left  by  we  obtain 


1  + 


m 


(B6) 


Similarly, 


1  +  C*R"^C 

m  o  a  c 


(B7) 


To  establish  the  boundedness  of  r(z)  it  is  necessary  to  prove  that 

Im-  r(z)r*{z)  >0  ,  |z|  <  1  .  (B8) 
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From  (B3),  (Bl),  and  (B2)  we  obtain,  ^ 

(1  -rr*)R"^l  -c  W’*‘)=(l  -WC*)R"^(1-C  W*)-(W-C  )Rr^(W*-C*) 

'm  o  a'm  a'mo''m  o'  a  o  o'  o'  o 

=  R"^-WRr^W*-C  R-'^C’*‘+WC*R"^C  W*  + 
a  D  ODO  oao 

+  C  Rr^W*-R^C  W*+WRr^C*-WC*R"^ 

OD  ao  DO  oa 

=  1  -WW*  .  (B9) 

m  '  ' 

Thus, 

1  -rr*  =  R(i  -wc*)"Vi  -ww*)(i  -c  w*)‘^R  >o  (bio) 

m  a'm  o''m  ''mo'  a  — 

since  1^-WW*>  0  (by  hypothesis),  Q.E.  D. 

The  proof  of  the  following  impedance  version  is  considerably  simpler 
and  is  therefore  omitted  [8]. 

Corollary.  Let  be  a  constant  mxm  matrix  such  that 

C  J- 

R  =— 2_°l  >  0  .  -;«ii) 

o  2 

l/2 

Let  R^'  equal  the  (unique)  hermitean  positive -definite  square-ixrot  of  R^. 

Then,  if  W(z)  is  an  mxm  positive  matrix, 

r(z)  =  R"^/^(W-C^)(W  +  C*)’^Ry^  (B12) 

is  an  mxm  bounded  matrix. 

NOTE:  By  inverting  (B12)  we  obtain 

c  •** 

W(z)  =  Ry^(l^+r)(l^-r)Ry^-l-  °g--  .  (B13) 

In  particular,  if  is  also  hermitean,  R^=  and 

w(z)  =  cy^i^+r)(i^-r)“^cy^  .  (B14) 

I  -2  -2 

By  starring  both  sides  of  (B5)  we  get  R^  C^=  C^R^  . 
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PART  IV 

Some  Observations  and 
Suggestions  for  New  Work 

by 

Haywood  E.  Webb  Jr. 


IV.  SOME  OBSERVATIONS  AND  SUGGESTIONS  FOR  NEW  WORK 


The  discussion  on  The  Interpolatory  Cascade,  in  Part  I  indicates 
that  every  admissible  spectral  estimator  can  be  interpreted  as  the  real 
part  of  the  driving  point  impedance  of  a  cascade  of  uniform  transmission 
lines  of  characteristic  impedances  given  by  (63)  terminated  in  an  arbitrary 
load  impedance  W(Z)  as  shown  in  Fig  1.  For  given  sample  covariance  data,  it 
is  a  different  load  impedance  that  corresponds  to  a  different  spectral 
estimator.  The  FEE  is  a  particular  class  of  estimators  which  provides  a 
set  of  "tuning"  parameters.  The  setting  of  all  the  work  here  is  in  a 
discrete  setting.  It  is  the  digital  computer  implementation  that  drives 
the  discrete  setting.  We  recall  that  a  process  with  a  continuous  covariance 
function,  in  order  to  be  non  determinstic,  must  satisfy  the  continuous 
analogue  of  (7)  Part  1. 

<  »  ^1) 

l+w'^ 


Where  A(wj  is  here  the  spectrum  obtained  from  the  Fourier  transform  of  the 
covariance.  But  note  also,  that  any  function  which  satisfies  (1)  is  not 
bandlimited,  so  that  a  discrete  sampling  representation  introduces  aliasing 
errors  at  the  outset,  and  the  problem,  from  a  continuous  viewpoint,  is 
improperly  formulated.  The  concept  of  entropy  does  not  make  sense. 

In  the  network  formulation  of  the  spectral  estimation  problem,  the 
discrete  cascade,  can  be  generalized  to  a  distributed  continuous  trans¬ 
mission  line.  Again,  we  conjecture,  that  every  admissible  spectrum 
obtainable  by  a  "continuation"  of  a  segment  of  the  covariance  function,  can 
be  obtained  by  a  non-uniform  transmission  line  with  a  shunt  capacitor, 
terminated  in  an  arbitrary  positive  real  load  impedance.  It  would  be 
interesting  to  maximize  the  continuous  analogue  (1)  in  this  section  to 
(7)  in  Part  1. 

Note  that  under  these  circumstances  that  (1)  of  Part  IV  cannot 
become  infinite  under  the  assumptions.  After  all,  fl)  of  Part  IV  is  a 
necessary  and  sufficient  condition  that  the  inverse  Fourier  transform  of 
.A(w)  not  be  quasi-analytic  i.e,,  not  be  uniquely  deterministic  over  its 
domain  from  its  derivatives  at  a  point. 

It  is  proposed  that  further  exploration  of  the  continuous  problem  be 
made,  and  that  after  understanding  it.  that  digital  computer  implementations 
be  considered  in  the  light  of  the  results. 
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